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Abstract. We study the astrophysical implications of neutralino dark matter annihilations in galaxy clusters, 
with a specific application to the Coma cluster. We first address the determination of the dark halo models for 
Coma, starting from structure formation models and observational data, and we discuss in detail the role of sub- 
halos. We then perform a thorough analysis of the transport and diffusion properties of neutralino annihilation 
products, and investigate the resulting multi-frequency signals, from radio to gamma-ray frequencies. We also 
study other relevant astrophysical effects of neutralino annihilations, like the DM-induced Sunyaev-Zel'dovich 
effect and the intracluster gas heating. As for the particle physics setup, we adopt a two-fold approach, resorting 
both to model-independent bottom-up scenarios and to benchmark, GUT-motivated frameworks. We show that 
the Coma radio-halo data (the spectrum and the surface brightness) can be nicely fitted by the neutralino-induced 
signal for peculiar particle physics models and for magnetic field values, which we outline in detail. Fitting the 
radio data and moving to higher frequencies, we find that the multi-frequency spectral energy distributions are 
typically dim at EUV and X-ray frequencies (with respect to the data), but show a non- negligible gamma-ray 
emission, depending on the amplitude of the Coma magnetic field. A simultaneous fit to the radio, EUV and HXR 
data is not possible without violating the gamma-ray EGRET upper limit. The best-fit particle physics models 
yields substantial heating of the intracluster gas, but not sufficient energy injection as to explain the quenching of 
cooling flows in the innermost region of clusters. Due to the specific multi-frequency features of the DM-induced 
spectral energy distribution in Coma, we find that supersymmetric models can be significantly and optimally 
constrained either in the gamma-rays or at radio and microwave frequencies. 
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1. Introduction 

Most of the matter content of the universe is in form of 
dark matter, whose presence is indicated by several as- 
trophysical evidences (e.g., gravitational lensing, galaxy 
rotation curves, galaxy clusters masses) but whose na- 
ture is still elusive. On the cosmological side, the most 
recent results of observational cosmology, i.e. WMAP vs. 
distant SN la, indicate that the matter content of the 
universe is Q m h 2 = 0.135lo'oo9 with a baryon density of 
n b h 2 = 0.0224 ± 0.0009 ( |Spergel et al. 2003| > . The combi- 
nation of the available data on large scale structures (Ly-a 
forest analysis of the SDSS, the SDSS galaxy clustering) 
with the latest SNe and with the 1-st year WMAP CMB 
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anisotropics can improve the determination of the cosmo- 
logical parameters ( |Seljak, et al. 2004| ) and hence allow us 
to set a concordance cosmological model. 
We refer, in this paper, to a flat ACDM cosmology with 
parameters chosen according to the global best fitting re- 
sults derived in [Seljak, et al. 2004| (see their Table 1, third 
column): we assume, in fact, that the present matter en- 
ergy density is VL m — 0.281, that the Hubble constant in 
units of 100 km s _1 Mpc -1 is h = 0.71, that the present 
mean energy density in baryons is fib = 0.0233//i 2 , with 
the only other significant extra matter term in cold dark 
matter flcDM — — ^b, that our Universe has a flat ge- 
ometry and a cosmological constant A, i.e. fl\ = 1 — fi m , 
and, finally, that the primordial power spectrum is scale 
invariant and is normalized to the value og = 0.897. This 



2 



S. Colafrancesco et al.: DM annihilations in Coma 



choice sets our framework, but it is not actually crucial 
for any of the results presented in the paper, which can 
be easily rescaled in case of a re- assessment of best-fit 
values of the cosmological parameters, and in particu- 
lar of the value of Sic dm (the present concordance cos- 
mological model, while widely used, has been also criti- 
cized and questioned in the light of still unexplored sys- 
tematics, see, e.g., |Myers et al. 2003) ISadat et al. 20051 
ILieu & Mittaz 20051 |Copi et al. 2003| ) 

The Coma cluster has been the first astrophysical lab- 
oratory for dark matter (DM) since the analysis of F. 
Zwicky ( |Zwicky T9 33). In this respect, we can consider the 
Coma cluster as an astrophysical benchmark case-study for 
DM. Modern observations have led to an increasingly so- 
phisticated exploration of the DM distribution in the uni- 
verse, now confirmed to be a dominant component (rel- 
ative to the baryonic material) over scales ranging from 
those of galaxy halos to that of the particle horizon. 

The nature of DM is not yet known and several detec- 
tion techniques have been used so far. Obviously, direct 
detection is the cleanest and most decisive discriminant 
(see e. a. IMunoz 20031 for a review). However, it would be 
interesting if astronomical techniques were to reveal some 
of the fundamental properties of DM particles. In fact, if 
DM is supposed to consist of fundamental particles which 
are weakly interacting, then their own interaction will lead 
to a number of astrophysical signatures (e.g., high-energy 
gamma-rays, electrons and positrons, protons and neutri- 
nos and hence by their emission/interaction properties) 
indicative of their nature and composition. 

These facts provide the basic motivations for our study, 
which is aimed to: i) describe the multi- wavelength signals 
of the presence of DM through the emission features of the 
secondary products of neutralino annihilation. These sig- 
nals are of non-thermal nature and cover the entire electro- 
magnetic spectrum, from radio to gamma-ray frequencies; 
ii) indicate the best frequency windows where it will be 
possible to catch physical indications for the nature of 
DM; in) apply this analysis to the largest bound contain- 
ers of DM in the universe, i.e. galaxy clusters. We shall 
focus here on the case of the Coma cluster, a particularly 
rich and suitable laboratory for which an extended obser- 
vational database is at hand. 

1.1. The fundamental physics framework 

Several candidates have been proposed as fundamental 
DM constituents, ranging from axions to light, MeV DM, 
from KK particles, branons, primordial BHs, mirror mat- 
ter to supersymmetric WIMPs (see, e.g., |Baltz, 2004| 
IBertone et al. 20041 and |Bergstrom 2000| for recent re- 
views). In this paper we will assume that the main DM 
constituent is the lightest neutralino of the minimal su- 
persymmetric extension of the Standard Model (MSSM). 
Although no experimental evidence in favor of supersym- 
metry has shown up to date, several theoretical moti- 
vations indicate that the MSSM is one of the best bets 



for new physics beyond the Standard Model. Intriguingly 
enough, and contrary to the majority of other parti- 
cle physics candidates for DM, supersymmetry can un- 
ambiguously manifest itself in future accelerator experi- 
ments. Furthermore, provided neutralinos are stable and 
are the lightest supersymmetric particles, next generation 
direct detection experiments feature good chances to ex- 
plore most of the neutralino-proton scattering cross sec- 
tion range predicted by supersymmetry. 

A long standing issue in phenomcnological studies 
of low-energy supersymmetry is traced to the param- 
eterization of the supersymmetry breaking terms (see 
|Chung et al. 200"3| for a recent review). In this respect, two 
somehow complementary attitudes have been pursued. On 
the one hand, one can appeal to a (set of) underlying high 
energy principles to constrain the form of supersymmetry 
breaking term, possibly at some high energy (often at a 
grand unification) scale (see e.g. IBaer et al. 2000|) . The 
low energy setup is then derived through the renormal- 
ization group evolution of the supersymmetry breaking 
parameters down to the electroweak scale. Alternatively, 
one can directly face the most general possible low energy 
realization of the MSSM, and try to figure out whether 
general properties of supersymmetry phenomenology can 
be derived (see e.g. |Profumo and Yaguna 2004| ). 

In this paper we will resort to both approaches. We 
will show that the final state products of neutralino pair 
annihilations show relatively few spectral patterns, and 
that any supersymmetric configuration can be thought as 
an interpolation among the extreme cases we shall con- 
sider here. The huge number of free parameters of the 
general MSSM are therefore effectively decoupled, and the 
only relevant physical properties are the final state prod- 
ucts of neutralino pair annihilations, and the mass of the 
neutralino itself. We will indicate this first strategy as a 
bottom-up approach (see Sect l4.l| for details). 

Since most phenomenological studies have been so far 
based on GUT-motivated models, and a wealth of re- 
sults on accelerator physics, direct and indirect detec- 
tion has accumulated within these frameworks, we de- 
cided to work out here, as well, the astrophysical con- 
sequences, for the system under consideration, of a few 
benchmark models. The latter have been chosen among 
the minimal supergravity (mSUGRA) models indicated 
in Batta gha et al. 2003| with the criterion of exemplifying 
the widest range of possibilities within that particular the- 
oretical setup (see Sect 13. 21 for details). 

1.2. The astrophysical framework 

To make our study quantitative, we will compare 
the predictions of the above mentioned neutralino 
models with the observational set-up of the Coma 
cluster, which represents the largest available observa- 
tional database for a galaxy cluster. The total mass 
of Coma found within 10ft. -1 Mpc from its center is 
M <wh -x Mpc « 1.65 x 1O 15 M feller et al. 1999|l. The 
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assumption of hydrostatic equilibrium of the thermal 
intra-cluster gas in Coma provides a complementary es- 
timate of its total mass enclosed in the radius r. A value 
M w 1.85 x 10 15 M Q within 5h~^ Mpc from the cluster 
center has been obtained from X-ray data ( |Hughes "19 89 ) . 
X-ray observations of Coma also yield detailed informa- 
tion about the thermal electrons population. We know 
that the hot thermal electrons are at a temperature 
k B T e = 8.2 ± 0.4 keV (|Arnaud et al. 2f)01|) and have 
a central density n = (3.42 ± 0.047) h 1 ^ 2 x 10" 3 
cm~ 3 , with a spatial distribution fitted by a f3- 
modcl, n(r) — no(l + r 2 /r 2 ) -3 ' 3 / 2 , with parameters 
r c = 10.5' ± 0.6' and /3 » 0.75 IjBriel et al. 1992|) . 
Assuming spherical symmetry and the previous parame- 
ter values, the optical depth of the thermal gas in Coma 
is T t h — 5.54 x 10~ 3 and the pressure due to the thermal 
electron population is P th ~ 2.80 • 10~ 2 keV cm" 3 . The 
hot intra-cluster gas produces also a thermal SZ effect 
HSunyaev fc Zel'dovich 19721 |Sunyaev fc ZeFdovich 1980| 
see|Birkinshaw 1999 for a general review) which has been 
observed over a wide frequency range, from 32 to 245 
GHz (see lDePetris et al. 20031 and references therein). 
Beyond the presence of DM and thermal gas, Coma 
also shows hints for the presence of relativistic particles 
in its atmosphere. The main evidence for the presence 
of a non-thermal population of relativistic electrons 
comes from the observation of the diffuse radio halo at 
frequencies v r ~ SOMHz — 5GHz (De iss et al. 19971 
IThierbach et al. "2 003 ) . The radio halo spectrum can 
be fitted by a power-law spectrum J v ~ j,— 1.35 m ^ ne 
range 30 MHz-1.4 GHz with a further steepening of the 
spectrum at higher radio frequencies. The radio halo of 
Coma has an extension of Rh ~ 0.9/ifQ 1 Mpc, and its 
surface brightness is quite flat in the inner 20 arcmin 
with a radial decline at larger angular distances (e.g., 



Other diffuse non-thermal emissions have been reported 
for Coma (as well as for a few other clusters) in the 
extreme UV (EUV) and in the hard X-ray (HXR) 
energy bands. The Coma flux observed in the 65 — 245 
eV band IjLieu et al. 19 96 ) is ~ 36% above the ex- 
pected flux from the thermal bremsstrahlung emission 
of the k B T w 8.2 keV IC gas l|Ensslin et al. 1998jl and 
it can be modeled with a power-law spectrum with 
an approximately constant slope ~ 1.75, in different 
spatial regions <|Lieu et al. 19991 lEnsslin et al. 1 999 
|Bowyer et al. 2004| ). The EUV excess in Coma has been 
unambiguously detected and it does not depend much on 
the data analysis procedure. The integrated flux in the 
energy band .13- 0.18 keV is F E UV « (4.1 ± 0.4) • 10~ 12 
erg cm _2 s _1 ( |Bowyer et al. 2004| |. According to the most 
recent analysis of the EUVE data ( |Bowyer et al. 200"4| ), 
the EUV excess seems to be spatially concentrated in 
the inner region (9 ^ 15 — 20 arcmin) of Coma (see 
also IBonamente et al. 2003JI . The nature of this excess 
is not definitely determined since both thermal and 
non-thermal models are able to reproduce the observed 



EUV flux. However, the analysis of |Bowyer et al. 2 004 
seems to favour a non-thermal origin of the EUV ex- 
cess in Coma generated by an additional population 
of secondary electrons. A soft X-ray (SXR) excess (in 
the energy range « 0.1 — 0.245 keV) has been also 
detected in the outer region (20' < 9 < 90') of Coma 
<|Bonamente et al. 20031 |Finoguenov et al. 2003) . The 
spectral features of this SXR excess seem to be more 
consistent with a thermal nature of the SXR emission, 
while a non-thermal model is not able to reproduce 
accurately the SXR data (e.g., IBonamente et al."2~0 03). 
The SXR emission from the outskirts of Coma has 
been fitted in a scenario in which the thermal gas at 
fcgT e ~ 0.2 keV with ~ 0.1 solar abundance (see, e.g. 
|Finoguenov et al. 2003| who identify the warm gas with a 
WHIM component) resides in the low-density filaments 
predicted to form around clusters as a result of the 
evolution of the large-scale web-like structure of the 
universe fsee lBonamente et al. 20031 . It nas been noticed, 
however, that the WHIM component cannot reproduce, 
by itself, the Coma SXR excess because it would produce 
a SXR emission by far lower (see IMittaz et al. 200H , and 
thus one is forced to assume a large amount of warm 
gas in the outskirts of Coma. Thus, it seems that the 
available EUV and SXR data indicate (at least) two 
different electron populations: a non-thermal one, likely 
yielding the centrally concentrated EUV excess and a 
thermal (likely warm) one, providing the peripherically 
located SXR excess. In this paper, we will compare our 
models with the EUV excess only, which is intimately 
related to Coma being spatially concentrated towards the 
inner region of the cluster. 

There is also evidence of a hard X-ray (HXR) emis- 
sion observed towards the direction of Coma with 
the BeppoSAX-PDS IjFusco-Femiano et al. 19991 
IFusco-Femiano et al. 2004| and w ith the ROSSI-XTE 
experiments (Rc phaeli et al. 1999| ). Both these measure- 
ments indicate an excess over the thermal emission which 
amounts to -F(2o-80)fc e y = ( 1-5 ± 0.5) ■ 10 -11 erg s _1 cm -2 
( Fus co-Femiano et al. 2004). It must be mentioned, for 
the sake of completeness, that the HXR excess of Coma 
is still controversial (see IRossetti and M olcndi 2004), 
but, for the aim of our discussion, it could at worst be 
regarded as an upper limit. The nature of the HXR 
emission of Coma is not yet fully understood. 
Finally, a gamma-ray upper limit of F(> lOOAfeU) « 
3.2 x 10~ 8 pho cm~ 2 has been derived for Coma 
from EGRET observations l|Sreekumar et al. 1 996 
IReimer et al. 2004jl . 

The current evidence for the radio-halo emission fea- 
tures of Coma has been interpreted as synchrotron emis- 
sion from a population of primary relativistic electrons 
which are subject to a continuous re-acceleration pro- 
cess supposedly triggered by merging shock and/or intra- 
cluster turbulence fe.g.. IBrunetti et al. 2004)l . The EUV 
and HXR emission excesses are currently interpreted 
as Inverse Compton scattering (ICS) emission from ei- 
ther primary or secondary electrons. Alternative mod- 
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eling has been proposed in terms of suprathermal elec- 
tron bremsstrahlung emission for the HXR emission of 
Coma (sec Ens slin et al. 19991 |Kempner fe Sarazin 2000| 
see also IPetrosian 20011 for a critical discussion) and in 
this case the EUV emission should be produced either by 
a different relativistic electron population or by a warm 
thermal population both concentrated towards the cluster 
center. Lastly, models in which the EUV and HXR emis- 
sion can be reproduced by synchrotron emission from the 
interaction of Ultra High Energy cosmic rays and/or pho- 
tons IjTimokhin et al. 20031 llnoue et al. 2005)l have also 
been presented. The situation is far from being completely 
clear and several problems still stand on both the obser- 
vational and theoretical sides of the issue. 

Since DM is abundantly present in Coma and relativis- 
tic particles are among the main annihilation products, we 
explore here the effect of DM annihilation on the multi- 
frequency spectral energy distribution (SED) of Coma. 
The plan of the paper is the following. We discuss in Sect.0 
the DM halo models for the Coma cluster, the set of best 
fitting parameters for the DM distribution, and the role of 
sub-halos. The annihilation features of neutralinos and the 
main annihilation products are discussed in Sect. |3| The 
multi-frequency signals of DM annihilation are presented 
and discussed in details in Sect. 0] and while the details 
of the transport and diffusion properties of the secondary 
particles are described in the Appendix lAl together with 
the derivation of the equilibrium spectrum of relativistic 
particles in Coma. The conclusions of our analysis and 
the outline for forthcoming astrophysical searches for DM 
signals in galaxy clusters are presented in the final Sect El 

2. A ACDM model for the Coma cluster 

To describe the DM halo profile of the Coma cluster we 
refer, as a general setup, to the ACDM model for struc- 
ture formation, implementing results of galaxy cluster for- 
mation obtained from N-body simulations. Free parame- 
ters are fitted against the available dynamical informa- 
tion and are compared to the predictions of this scheme. 
Substructures will play a major role when we will discuss 
the predictions for signals of DM annihilations. In this re- 
spect, the picture derived from simulations is less clean 
and, hence, we will describe in details our set of assump- 
tions. 



2.1. The dark matter halo profile 

To describe the DM halo profile of the Coma cluster we 
consider the limit in which the mean DM distribution in 
Coma can be regarded as spherically symmetric and rep- 
resented by the parametric radial density profile: 



P(r) = p'g(r/a) 



(1) 



Two schemes are adopted to choose the function g{x) in- 
troduced here. In the first one, we assume that g(x) can 
be directly inferred as the function setting the universal 



shape of DM halos found in numerical N-body simula- 
tions of hierarchical clustering. We are assuming, hence, 
that the DM profile is essentially unaltered from the stage 
preceding the baryon collapse, which is - strictly speaking 
- the picture provided by the simulations for the present- 
day cluster morphology. A few forms for the universal DM 
profile have been proposed in the literature: we implement 
here the non-singular form (which we label as N04 profile) 
extrapolated bv lNavarro et al.~2 004 



<7Aro4(a:) = exp[— 2/a{x a — 1)] with 



0.17 , 



(2) 



and the shape with a mild singularity towards its cen- 
ter proposed bv lDiemand et al. 20051 (labeled here as D05 
profile): 



9do5(x) 



\3- 7 



with 7 ~ 1.2 



(3) 



The other extreme scheme is a picture in which the baryon 
infall induces a large transfer of angular momentum be- 
tween the luminous and the dark components of the cos- 
mic structure, with significant modification of the shape 
of the DM profile in its inner region. According to a re- 
cent model IjEl-Zant et al. 2001|l . baryons might sink in 
the central part of DM halos after getting clumped into 
dense gas clouds, with the halo density profile in the fi- 
nal configuration found to be described by a profile (la- 
beled here as B profile) with a large core radius (see, e.g., 
IBurkert 1995) : 



9b(x) 



(4) 



Once the shape of the DM profile is chosen, the radial 
density profile in Eq. Q is fully specified by two param- 
eters: the length-scale a and the normalization parameter 
p' . It is, however, useful to describe the density profile 
model by other two parameters, i.e., its virial mass M V i r 
and concentration parameter c V i r . For the latter param- 
eter, we adopt here the definition bv lBullock et al. 20011 
We introduce the virial radius R V i r of a halo of mass M V i r 
as the radius within which the mean density of the halo is 
equal to the virial overdensity times the mean back- 
ground density p = fl m pc : 



_ 4tt 



pR 



3 

vir ' 



(5) 



We assume here that the virial overdensity can be approx- 
imated by the expression (see |Bryan fc Norman~9 98 ) , 
appropriate for a flat cosmology, 



(18tt 2 + 82x- 39x 2 ) 
1-x 



(6) 



with x = Vt m (z) — 1. In our cosmological setup we find at 
z = 0, A vir ~ 343 (we refer to IColafrancesco et al. 1 994 
ICola fraiicesc o et al. 19971 for a general derivation of the 
virial overdensity in different cosmological models). The 
concentration parameter is then defined as 



R,, 



R, 



r-2 



X-2 a 



(7) 
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with r_2 the radius at which the effective logarithmic slope 
of the profile is -2. We find that x_ 2 = 1 for the N04 
profile (see Eq.[2Jl, x^ 2 — 2 — 7 for D05 profile (see Eq.|3), 
and x_2 — 1.52 for the Burkert profile (see Eq. 0}. 

Since the first numerical results with large statistics 
became available IjNavarro et al. 19 97). it has been real- 
ized that, at any given redshift, there is a strong corre- 
lation between c V i r and M V i r , with larger concentrations 
found in lighter halos. This trend may be intuitively ex- 
plained by the fact that mean overdensities in halos should 
be correlated with the mean background densities at the 
time of collapse, and in the hierarchical structure forma- 
tion model small objects form first, when the Universe was 
indeed denser. The correlation between c V i r and M V i r is 
relevant in our context at two levels: i) when discussing 
the mean density profile of Coma and, ii) when including 
substructures. Hence, we will review this relevant issue 
here and we will apply it to the present case of Coma. 
Bull ock et al. 200"T1 proposed a model to describe this cor- 
relation, improving on the toy model originally outlined 
in INavarro et al. 19971 A collapse redshift z c is assigned, 
on average, to each halo of mass M at the epoch z through 
the relation M*(z c ) = FM. Here it is postulated that a 
fixed fraction F of M (following Wcch sler et al. 20011 we 
choose F = 0.015) corresponds to the typical collapsing 
mass M*, as defined implicitly by a (M±(z)) — 5 sc (z), with 
5 SC being the critical overdensity required for collapse in 
the spherical model and a(M) being the present-day rms 
density fluctuation in spheres containing a mean mass M 
(see, e.g., IPeebles 1 980). An expression for 5 SC is given, 
e.g., in lEke et al. 19961 The rms fluctuation a{M) is re- 
lated to the fluctuation power spectrum P(k) (see e.g. 
IPeehles 1993) by 



a 2 {M) = d 3 kW 2 (kR)P{k). 



(8) 



where W is the top-hat window function on the scale 
R 3 = 3M/4np with p the mean (proper) matter density, 
i.e. p = £l m Pc with p c the critical density. The power 
spectrum P(k) is parametrized as P(k) oc k n T 2 (k) in 
terms of the primordial power-spectrum shape oc k n 
and of the transfer function T 2 (k) associated to the 
specific DM scenario. We fix the primordial spectral 
index n = 1 and we take the transfer function T 2 (k) 
given by IBardeen et al. 1 986 for an adiabatic CDM 
model, with the shape parameter modified to include 
baryonic matter according to the prescription in, 
e.g. IPeacock 19991 (see their eqs. 15.84 and 15.85) and 
introducing a multiplicative exponential cutoff at large 
k corresponding to the free-streaming scale for WIMPs 
Pofmann et al. 2001, Chen et al. 2001, Creen et al. 2005~ 
The spectrum P(k) is normalized to the value as = 0.897 
as was quoted above. 

The toy model of Bullock et al. (2001) prescribes a one 
to one correspondence between the density of the Universe 
at the collapse redshift z c of the DM halo and a character- 
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Fig. 1. The dependence of c V i r on the halo mass M, at z = 
0, as in the Bullock et al. toy model (solid line) and in the 
ENS toy model (dashed line) ; predictions are compared to 
a few sets of simulation results in different mass ranges. A 
flat, vacuum-dominated cosmology with VLm — 0.3, Qa = 
0.7, h — 0.7 and as = 1 is assumed here. 

istic density of the halo at the redshift z; it follows that, 
on average, the concentration parameter is given by 



,{M,z) 



K 



r l + z c _ c mr (M, z = 0) 



1 



1 + 2) 



(9) 



with K being a constant (i.e. independent of M and 
cosmology) to be fitted to the results of the N-body 
simulations. We plot in Fig. the dependence of 
c v i r on the halo mass M, at z = 0, according to 
the toy model of Bullock et al. (2001) as extrap- 
olated down to the free-streaming mass scale for 
DM halos made of WIMPs, i.e. around 10" 6 M Q (see 
Hof mann et al. 2001, Chen et al. 2001, Green et al. 2005, Diemand et 
The predictions are compared to the results of a few 
sets of N-body simulations: we use "data" points and 
relative error bars from IBullock et al. 200D (represent- 
ing a binning in mass of results for a large sample of 
simulated halos; in each mass bin, the marker and the 
error bars correspond, respectively, to the peak and the 
68% width in the c V i r distribution) to determine the 
parameter K. The same value will be used to infer the 
mean c V i r predicted in our cosmological setup. Other 
"datasets" refer actually to different values of as and 
diff erent redshifts z (z = 26 for the two minihalos fitted 
r^gW^^Sic^HM ct a i. 2005 and for the upper bound 
in the range up to 10 Mq quoted in the same paper; 
z = 3 for the sample from IColin et al. 2004]! and have 
been extrapolated, consistently with our prescriptions, to 
z = and erg = 1. Since small objects tend to collapse 
all at the same redshift, the dependence on mass of the 
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concentration parameters flattens at small masses; the 
mean asymptotic value we find is slightly larger than the 
typical values found in D iemand et al. 20051 but it is still 
consistent with that analysis. 

An alternative toy-model to describe the relation be- 
tween c v i r and M has been discussed by Eke, Navarro and 
Steinmetz IjEke et al. 20011 hereafter ENS model). The re- 
lation they propose has a similar scaling in z, but with a 
different definition of the collapse redshift z c and a milder 
dependence of c V i r on M. In our notation, they define z c 
through the equation 



D(z c )a eS (M p ) = — 



(10) 



where D(z) represents the linear theory growth factor, and 
cr c ff is an 'effective' amplitude of the power spectrum on 
scale M: 



a eS (M) = a(M) 



dln(cr) 



(M) 



da 
dM' 



-M 



(11) 



d\n(M) 

which modulates <r(M) and makes z c dependent on both 
the amplitude and on the shape of the power spectrum, 
rather than just on the amplitude, as in the toy model 
of Bullock et al. (2001). Finally, in Eq. JHI, M p is as- 
sumed to be the mass of the halo contained within the 
radius at which the circular velocity reaches its maxi- 
mum, while C<j is a free parameter (independent of M 
and cosmology) which we will fit again to the "data" set 
in Bullo ck~et al. 200T1 With such a definition of z c it fol- 
lows that, on average, c V i r can be expressed as: 



•(M,z) 



A vir (z c ) D M (z) 



1/3 



1 + Zc 



(12) 



A V i r (z)Q M (zc) J l + z 

As shown in Fig. ^ the dependence of c„j r on M given 
by Eci. (|12[) above is weaker than that obtained in the 
Bullock et al. (2001) toy-model, with a significant mis- 
match in the extrapolation already with respect to the 
sample from IColin et al. 20041 and an even larger mis- 
match in the low mass end. Moreover, the extrapolation 
breaks down when the logarithmic derivative of the a(M) 
becomes very small, in the regime when P{k) scales as 
fc~ 3 . Note also that predictions in this model are rather 
sensitive to the specific spectrum P(k) assumed (in par- 
ticular the form in the public release of the ENS numerical 
code gives slightly larger values of c V i r in its low mass end, 
around a value c V i r « 40 (we checked that implementing 
our fitting function for the power spectrum, we recover 
our trend). 

2.2. Fitting the halo parameters of Coma 

For a given shape of the halo profile we make a fit of the 
parameters M V i r and c V i r against the available dynamical 
constraints for Coma. We consider two bounds on the to- 
tal mass of the cluster at large radii, as inferred with tech- 
niques largely insensitive to the details of the mass profile 
in its inner region. In Gell er et al. 19 99. a total mass 

M(r < 10 IT 1 Mpc) = (1.65 ± 0.41) 10 15 h" 1 M (13) 



is derived mapping the caustics in redshift space of galax- 
ies infalling in Coma on nearly radial orbits. Several 
authors derived mass budgets for Coma using optical 
data and applying the virial theorem, or using X-ray 
data and assuming hydrostatic equilibrium. We consider 
the bound derived by |Hughes 19 89, cross-correlating such 
techniques: 



(14) 



M(r < 5 h^ Mpc) = (1.85 ± 0.25) 10 15 h^ 1 M 

where h^o is the Hubble constant in units of 50 km s _1 
Mpc" 1 . 

In our discussion some information on the inner shape 
of the mass profile in Coma is also important: we imple- 
ment here the constraint that can derived by studying 
the velocity moments of a given tracer population in the 
cluster. As the most reliable observable quantity one can 
consider the projection along the line of sight of the radial 
velocity dispersion of the population; under the assump- 
tion of spherical symmetry and without bulk rotation, this 
is related to the total mass profile M(r) by the expression 
UBinney fe Mamon 1 982"! ILokas & Mamon 2003|> : 

<W = £Lr dr ' v{ri)M{r'){r>r- 2 



i(R) Jr 



-2/3+1 



vr^T 



R 2 ' 



(15) 



where v(r) is the density profile of the tracer population 
and I(R) represents its surface density at the projected 
radius R. In the derivation of Eq. I|15|l . a constant-over- 
radius anisotropy parameter (3 defined as 



(3 = 1- 



a 2 (r) 



(16) 



has been assumed with a 2 and cr| being, respectively, the 
radial and tangential velocity dispersion {[3 = 1 denotes 
the case of purely radial orbits, (3 = that of system 
with isotropic velocity dispersion, while (3 — > — oo labels 
circular orbits). Following ILokas fc Mamon 20031 we take 
as tracer population that of the E-S0 galaxies, whose line 
of sight velocity dispersion has been mapped, according 
to Gaussian distribution, in nine radial bins from 4' out 
to 190' (see Fig. 3 in ILokas fc Mamon 2003|l . and whose 
density profile can be described by the fitting function: 



1 



j/ r s)(l + r/r s ) 2 



(17) 



with ?'s = 7 .'05. Constraints to the DM profile are obtained 
through its contribution to M(r), in which we include the 
terms due to spiral and E-S0 galaxies (each one with the 
appropriate density profile normalized to the observed lu- 
minosity through an appropriate mass-to-light ratio), and 
the gas component (as inferred from the X-ray surface 
brightness distribution) whose number density profile can 
be described by the fitting function: 



n(r) = no 



1+ - 



-1.5 b 



(18) 
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, r c = W.5 and b = 0.75 



with n = 3.42 x 1(T 3 cm" 3 
IjBriel et al. 1992jl . 

To compare a model with such datasets, we build a 
reduced x 2 -like variable of the form: 



X r 



2 (M(r<n)-Mi) 2 



(AM, 



i=l 

9 ^ 



(A< s ) 2 



(19) 



where the index i in the first sum runs over the constraints 
given in Eqs. 1|13|1 and l|14|l. while, in the second sum, we 
include the nine radial bins over which the line of sight 
velocity dispersion of E-SO galaxies and its standard devi- 
ation has been estimated. Weight factors have been added 
to give the same statistical weight to each of the two 
classes of constraints, see, e.g., |Dehnen fc Binne~1 998 
where an analogous procedure has been adopted. 

Nonetheless, we have derived in Fig. |3 the 1 a, 2 a 
and 3 a contours in the (M V i r , c V i r ) plane for the Navarro 
et al. halo profile (Eq. © and for the Diemand et al. 
halo profile (Eq.|2J). In Fig.|31we show the analogous con- 
tours for the Burkert profile (Eq. JJJl). In all these cases we 
have performed the fit of the line-of-sight radial velocity 
dispersion of E-SO galaxies assuming that this system has 
an isotropic velocity dispersion, i.e. we have taken (3 = 0. 
Best fitting values are found at M v i r ~ 0.9 • 1O 15 M hr 1 
and c v i r ~ 10 (that we consider, hence, as reference values 
in the following analysis), not too far from the mean value 
expected from models sketching the correlation between 
these two parameters in the ACDM picture. We show in 
Figs. [5] and [3] the predictions of such correlation in the 
models of Bullock et al. (solid line) and of lEke et al. 2 001 
(dashed line). 

2.3. Substructures in the Coma cluster 

Since the astrophysical signals produced by WIMP pair 
annihilation scale with the square of the WIMP den- 
sity, any local overdensity does play a role (see e.g. 
|Bergstrom et al. 1998| and references therein). To discuss 
substructures in the Coma cluster, analogously to the gen- 
eral picture introduced above for DM halos, we label a 
subhalo through its virial mass M s and its concentration 
parameter c s (or equivalently a typical density and length 
scale, p' s and a s ). The subhalo profile shape is considered 
here to be spherical and of the same form as for the par- 
ent halo. Finally, as for the mean DM density profile, the 
distribution of subhalos in Coma is taken to be spherically 
symmetric. The subhalo number density probability dis- 
tribution can then be fully specified through M s , c s and 
the radial coordinate for the subhalo position r. To our 
purposes, it is sufficient to consider the simplified case 
when the dependence on these three parameters can be 
factorized, i.e.: 



dr 3 dM„ dc. 



(20) 
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Fig. 2. We show the 1 cr, 2 a and 3 a contours as derived 
from the Xr variable in eq.(19), for the Navarro et al. halo 
profile, Eq. (J2J, (upper panel) and for the Diemand et al. 
halo profile, Eq. (0 (lower panel). Also shown are mean 
values for the correlation between M V i r and c V i r as in the 
toy models of Bulloc k~et al. 200~T1 (solid line) and that of 
lEke et al. 200D Mashed line). 



Here we have introduced a subhalo mass function, inde- 
pendent of radius, which is assumed to be of the form: 



dn s 
dM~ s 



A(M V 



Mi 



exp 



M cu 



-2/3' 



(21) 



IDiemand et al. 20051 

where M cut is the free streaming cutoff mass 

plofmann et al. 2001, Chen et al. 2001, Creen et al. 2005, Diemand et 

while the normalization A(M V i r ) is derived imposing that 
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0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 

M vjr [10 15 M h" 1 ] 

Fig. 3. We show here the 1 a, 2 a and 3 a contours as 
derived from the xl variable introduced in the text, for 
the Burkert profile, Eq. (gj. 



the total mass in subhalos is a fraction f s of the total 
virial mass M vir of the parent halo, i.e. 



M vir , 

dMs ^M s 
s dM s 



fsMv 



(22) 



According to IDiemand et al. 20051 f s is about 50% for a 
Milky Way size halo, and we will assume that the same 
holds for Coma. The quantity Vsics) is a log-normal distri- 
bution in concentration parameters around a mean value 
set by the substructure mass; the trend linking the mean 
c s to M s is expected to be analogous to that sketched 
above for parent halos with the Bullock et al. or ENS toy 
models, except that, on average, substructures collapsed 
in higher density environments and suffered tidal strip- 
ping. Both of these effects go in the direction of driving 
larger concentrations, as observed in the numerical sim- 
ulation of Bulloc k~et al. 200~T1 where it is shown that, on 
average and for M ~ 5 • 10 n Mo objects, the concentra- 
tion parameter in subhalos is found to be a factor of rs 1.5 
larger than for halos. We make here the simplified ansazt: 



(c s (AT, )) = F s (c mr (M mr )) with M s = M v 



(23) 



where, for simplicity, we assume that the enhancement 
factor F s does not depend on M s . Following again 
Bull ock et al. 200ll the la deviation A(log 10 c s ) around 
the mean in the log- normal distribution V s (c s ), is assumed 
to be independent of M s and of cosmology, and to be, nu- 
merically, A(log 10 c s ) = 0.14. 

Finally, we have to specify the spatial distribution of sub- 
structures within the cluster. Numerical simulations, trac- 
ing tidal stripping, find radial distributions which are sig- 
nificantly less concentrated than that of the smooth DM 



component. This radial bias is introduced here assuming 
that: 

p s (r) oc g(r/a') , (24) 

with g being the same functional form introduced 
above for the parent halo, but with a' much larger 
than the length scale a found for Coma. Following 
Nag ai fc Kravtsov 20051 we nx a '/ a — Since the frac- 
tion / s of DM in subhalos refers to structures within the 
virial radius, the normalization of p s (r) follows from the 
requirement 

f-Rvir 

-in / r 2 Ps (r) = 1 . (25) 



3. Neutralino annihilations in Coma 

3.1. Statistical properties 

Having set the reference particle physics framework and 
specified the distribution of DM particles, we can now in- 
troduce the source function from neutralino pair annihila- 
tions. For any stable particle species i, generated promptly 
in the annihilation or produced in the decay and fragmen- 
tation processes of the annihilation yields, the source func- 
tion Qi(r,E) gives the number of particles per unit time, 
energy and volume element produced locally in space: 



Q l {r, E) 



f 



dN[ 
dE 



(E)B f M pails (r) , 



(26) 



where (crw)o is the neutralino annihilation rate at zero 
temperature. The sum is over all kinematically allowed 
annihilation final states /, each with a branching ratio 
Bf and a spectral distribution dN? /dE, and A/" pa i rs (r) is 
the number density of neutralino pairs at a given radius 
r (i.e., the number of DM particles pairs per volume el- 
ement squared). The particle physics framework sets the 
quantity (erv)o and the list of Bf. Since the neutralino 
is a Majorana fermion light fermion final states are sup- 
pressed, while - depending on mass and composition - the 
dominant channels are either those with heavy fermions 
or those with gauge and Higgs bosons. The spectral func- 
tions d/v/ / dE are inferred from the results of MonteCarlo 
codes, namely the Pythia dSjostrand 1994, 1995[) 6.154, a s 
included in the DarkSUSY package HGondolo et al. 2004(1 . 
Finally, Af pa i rs (r) is obtained by summing the contribution 
from the smooth DM component, which we write here as 
the difference between the cumulative profile and the term 
that at a given radius is bound in subhalos, and the con- 
tributions from each subhalo, in the limit of unresolved 
substructures and in view of the fact that we consider 
only spherically averaged observables: 

Kr /s ~{p'g{r/a)~ f s M mrPs (r)f 
■A/pairs \T) = 



+ Ps (r) / dM : 



2 Ml 
dn„ 



'dM, 



dc' s V s (c s '(AQ) 



d 3 i 



{p' s 9( r s/a s )Y 
2AQ 



(27) 
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This quantity can be rewritten in the more compact form: 



A/pairs (?") — 



2M 2 



-f s A- 



(p'g(r/a) - f s p s g{r/a')) 



Psg(r/a') 



(28) 



where we have normalized densities to the present-day 
mean matter density in the Universe p, and we have de- 
fined the quantity: 



/.A 



with 



2 _ 



fdM s ^-M s A 2 M (M s ) 



fdM s -^-M s A 2 Ms (M s ) 



A 2 Ms (Ms) 



and 



I n {x) 



-(*) 



dn s 
<1M.< 



dc' s V s (c' B ) 



I 2 (c>_ 



\h{c' s x- 2 )] 



(29) 
(30) 

:{c' s X- 2 f (31) 




10- 



10 7 10' 



10 3 10" 



10 



io 3 io s 



7 9 1 

10 10 10 



io 13 io 15 



M [M, 



dyy 2 [g(y)] r 



(32) 



Such definitions are useful since A 2 M gives the average 
enhancement in the source due to a subhalo of mass M s , 
while A 2 is the sum over all such contributions weighted 
over the subhalo mass function times mass. Finally, in 
Eq. (|28|l we have also introduced the quantity: 



M„ 



4i(a') 3 Ji(E„/a') 



(33) 



In the limit in which the radial distribution of substruc- 
tures traces the DM profile, i.e. a' — a, p s becomes equal 
to the halo normalization parameter p' . 

We show in Fig. 0] the scaling of the average enhance- 
ment A 2 f in the source function versus the subhalo mass 
M s . We have considered the three halo models introduced 
in the previous Section, i.e. the N04, D05 and Burkert 
profiles, for the two toy models describing the scaling of 
concentration parameter with mass, i.e. the Bullock et al. 
and the ENS schemes, as well as two sample values for 
the ratio F s between the average concentration parame- 
ter in subhalos and that in halos of equal mass. In each 
setup, going to smaller and smaller values of M Sl the aver- 
age enhancement A\j increases and then flattens out at 
the mass scale below which all structures tend to collapse 
at the same epoch, and hence have equal concentration 
parameter. 

In Fig. [3] we show the scaling of the weighted enhance- 
ment A 2 in the source function due to subhalos versus 
the ratio between concentration parameter in subhalos to 
concentration parameter in halos at equal mass F s ; we 
give results for the usual set of halo profiles considered in 
our approach. Analogously to the enhancement for a fixed 
mass shown in the previous plot, A 2 is very sensitive to 
the scaling of the concentration parameter and hence we 



Fig. 4. Scaling of the average enhancement in source func- 
tions due to a subhalo of mass M s . We show result imple- 
menting the three halo profiles introduced, i.e. the N04, 
D05 and Burkert profile, the two toy models for the scal- 
ing of concentration parameter with mass, i.e. the Bullock 
et al. and the ENS, and two sample values of the ratio 
between concentration parameter in subhalos to concen- 
tration parameter in halos at equal mass F s . 



find a sharp dependence of A 2 on F s . The fractional con- 
tribution per logarithmic interval in subhalo mass M s to 
A 2 is also shown in Fig.[5]for four sample cases. Note that, 
although the factorization in the probability distribution 
for clumps in the radial coordinate and mass (plus the as- 
sumption that F s does not depend on mass) are a crude 
approximation, what we actually need in our discussion is 
F s and the radial distribution for subhalos at the peak of 
the distribution shown in Fig.|SJ unfortunately we cannot 
read out this from numerical simulations. 

Fig. [5] shows the number density of neutralino pairs 
(we set here the neutralino mass to M x = 100 GeV) as 
a function of the distance from the center of Coma for 
the three representative halo profiles introduced here, i.e. 
the N04, D05 and Burkert profile in their best fit model, 
and a sample configuration for the subhalo parameters. 
For the D05 and N04 profiles, the central enhancement 
increases the integrated source function by a factor « 6 
with respect to the Burkert profile, but this takes place 
on such a small angular scale that from the observational 
point of view it is like adding a point source at the center 
of the cluster. The enhancement of the annihilation signals 
from subhalos comes instead from large radii. This means 
that the enhancement from subhalos largely influences the 
results when the neutralino source is extended. This is the 
case of galaxy clusters, and more specifically of the Coma 
cluster which is our target in this paper. 
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Fig. 5. Left: scaling of the weighted enhancement in source functions due to subhalos versus the ratio between con- 
centration parameter in subhalos to concentration parameter in halos at equal mass F s . Results are shown the three 
halo profiles introduced, i.e. the N04, D05 and Burkert profile, the two toy models for the scaling of concentration 
parameter with mass, i.e. the Bullock et al. and the ENS. Right: fractional contribution per logarithmic interval in 
subhalo mass M s to A 2 in four sample cases. A normalization of the fluctuation spectrum tr 8 = 0.89 is adopted. 
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Fig. 6. The number density of neutralino pairs (neutralino 
mass set to M x — 100 GeV) as a function distance for the 
center of Coma for the three halo profiles introduced, i. e. 
the N04, D05 and Burkert profile in their best fit model, 
and a sample configuration for the subhalo parameters. 

3.2. Source functions spectral properties: generalities 
and supersymmetric benchmarks 

The spectral properties of secondary products of DM anni- 
hilations depend only, prior to diffusion and energy losses, 
on the DM particle mass M x and on the branching ra- 
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Fig. 7. The spectral shape of the electron source function 
in case of three sample final states (see text for details. 

tio BR(xx ~ ¥ f) f° r the final state / in the DM pair- 
annihilation. The DM particle physics model further sets 
the magnitude of the thermally averaged pair annihila- 
tion cross section times the relative DM particles velocity, 
(crv)o at T = 0. 

The range of neutralino masses and pair annihilation 
cross sections in the most general supersymmetric DM 
setup is extremely wide. Neutralinos as light as few GeV 
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(see lBottino et al. 2003(1 and as heavy as hundreds of TeV 
(see IProfumo 20 05 ) can account for the observed CDM 
density through thermal production mechanisms, and es- 
sentially no constraints apply in the case of non-thermally 
produced neutralinos. 

Turning to the viable range of neutralino pair an- 
nihilation cross sections, coannihilation processes do 
not allow us to set any lower bound, while on 
purely theoretical grounds a general upper limit on 

(cv)o ~ 10~ 22 (t^v) cm 3 /s has been recently set 
UProfumo 2005|l . The only general argument which ties the 
relic abundance of a WIMP with its pair annihilation cross 
section is given by the naive relation 



2 _ 3 x 10" 27 cm 3 /s 
(o-v)o 



(34) 



(sec Jung man et al. 1995| Eq.3.4), which points at a fidu- 
cial value for (av)o ~ 3 x 10 _26 cm 3 /s for our choice of cos- 
mological parameters. The above mentioned relation can 
be, however, badly violated in the general MSSM, or even 
within minimal setups, such as the minimal supergravity 
scenario ('see IProfumo 2005f . 

Since third generation leptons and quarks Yukawa cou- 
plings are always much larger than those of the first two 
generations, and being the neutralino a Majorana fermion, 
the largest BR(xx ~^ f) f° r annihilations into a fermion- 
antifcrmion pair are in most cases 1 into the third genera- 
tion final states bb, tt and t + t~. In the context of super- 
symmetry, if the supersymmetric partners of the above 
mentioned fermions are not significantly different in mass, 
the t + t~ branching ratio will be suppressed, with respect 
to the bb branching ratio by a color factor equal to 1/3, 
plus a possible further Yukawa coupling suppression, since 
the two final states share the same SU(2) quantum num- 
ber assignment. Further, the fragmentation functions of 
third generation quarks are very similar, and give rise to 
what we will dub in the following as a 11 soft spectrum!'' . 
A second possibility, when kinematically allowed, is the 
pair annihilation into massive gauge bosons 2 , W + W~ and 
Z{)Z{). Again, the fragmentation functions for these two 
final states are mostly indistinguishable, and will be in- 
dicated as giving a "hard spectrum" . The occurrence of a 
non-negligible branching fraction into t + t~ or into light 
quarks will generally give raise to intermediate spectra be- 
tween the "hard" and "soft" case. 

Figdshows the spectral shape of the electron source func- 
tion in the case of the three sample final states bb, t + t~ 
and W + W~ for M x — 100 GeV, and clarifies the previous 
discussion. In what follows we will therefore employ sam- 
ple DM configurations making use of cither soft (bb) or 

1 Models with non-universal Higgs masses at the GUT scale 
can give instances of exceptions to this generic spectral pattern, 
featuring light first and second generation sfermions (see e.g. 
IBaer et al. 2005bl . 

2 The direct annihilation into photons is loop suppressed in 
supersymmetric models (see e.g. |Bergstrom & Snellman 1988 
and |Bergstrom fc Ullio 1997| . 



hard (W^ + W ) spectra, keeping in mind that other possi- 
bilities would likely fall in between these two extrema. 

In order to make a more stringent contact with super- 
symmetry phenomenology, we will however also resort to 
realistic benchmark SUSY models: by this we mean thor- 
oughly defined SUSY setups which are fully consistent 
with accelerator and other phenomenological constraints, 
and which give a neutralino thermal relic abundance ex- 
actly matching the central cosmologically observed value. 
To this extent, we refer to the so-called minimal su- 
pergravity model HGoldberg 1983; Ellis et al. 1983, 1984| ), 
perhaps one of the better studied paradigms of low- 
energy supersymmetry, which enables, moreover, a cross- 
comparison with numerous dedicated studies, rang- 
ing from colliders (Bae r et al. 2003|) to DM searches 
jEdsjo et al. 2004, Baer et al. 20041 ). 

The assumptions of universality in the gaugino and 
in the scalar (masses and trilinear couplings) sectors re- 
markably reduce, in this model, the number of free pa- 
rameters of the general soft SUSY breaking Lagrangian 
| |Chung et al. 20 03) down to four continuous parame- 
ters (mo, Mi/ 2 , Ao, tan (3) plus one sign (sign(^)). The 
mSUGRA parameter space producing a sufficiently low 
thermal neutralino relic abundance f2 x /i 2 has been shown 
to be constrained to a handful of "regions" featuring ef- 
fective Q x h 2 suppression mechanisms (Ell is et al. 20 03). 
The latter are coannihilations of the neutralino with the 
next-to-lightest SUSY particle ( "Coannihilation" region) , 
rapid annihilations through s channel Higgs exchanges 
("Funnel" region), the occurrence of light enough neu- 
tralino and sfermions masses ("Bulk" region) and the 
presence of a non-negligible bino-higgsino mixing ( "Focus 
Point" region). 

With the idea of allowing a direct comparison with 
the existing research work in a wealth of complementary 
fields, we restrict ourselves to the "updated post- WMAP 
benchmarks for supersymmetr^ proposed and studied by 
Battagli a et al. 2003| All of those setups are tuned so as 
to feature a neutralino thermal relic density giving exactly 
the central WMAP-estimated CDM density 3 . As a prelim- 
inary step, we computed the electrons, neutrinos, gamma- 
rays and protons source spectra for all the 13 A'-M' mod- 
els. Remarkably enough, although the SUSY particle spec- 
trum is rather homogeneous throughout the mSUGRA pa- 
rameter space, the resulting spectra exhibit at least three 
qualitatively different shapes, according to the dominant 
final state in neutralino pair annihilation processes. In par- 
ticular, in the Bulk and Funnel regions the dominant final 
state is into bb, and, with a sub-dominant variable con- 
tribution, t + t~. The latter channel is instead dominant, 
for kinematic reasons, in the stau Coannihilation region. 
Finally, a third, and last, possibility is a dominant gauge 
bosons final state, which is the case along the Focus Point 
region. In this respect, in the effort to reproduce all of 
the mentioned spectral modes, and to reflect every cosmo- 



3 We adjusted here the values of mo given in 
Battaglia et al. 2003| in order to fulfill this requirement 
making use of the latest Isajet v. 7. 72 release and of the 
DarkSUSYpackage | |Edsjo et al. 2003| see Table[TJ. 
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Fig. 8. (Left): The electrons flux (dN e /dE)(av)o as a function of the electron energy. (Right): The gamma-rays flux 
(dN- f /dE)(av)o as a function of the photon energy. 




Fig. 9. (Left): The neutrinos flux (dN v / 'dE)(av)o as a function of the neutrino energy. (Right): The protons flux 
(dN p /dE)(av)o as a function of the proton energy. 



logically viable mSUGRA region, we focused on the four 
mod els indicated in Tab le ^ a subset of the benchmarks 
of Batt agha et al. 2003| (to which we refer the reader for 
further details). We collect in Table [5]the branching ratios 
for the final states of neutralino pair annihilations. In the last 
column of this table we also provide the thermally- averaged 
pair annihilation cross section times the relative velocity, at 
T — 0, {crv)o. Table|5|is an accurate guideline to interpret the 
resulting source spectra for the four benchmarks under consid- 



Model 


Ml/2 


mo 


tan/3 


sign(/x) 


rot 


B' (Bulk) 


250 


57 


10 


> 


175 


D' (Coann.) 


525 


101 


10 


> 


175 


E' (Focus P.) 


300 


1653 


10 


> 


171 


K' (Funnel) 


1300 


1070 


46 


< 


175 



Table 1. The input parameters of the four mSUGRA 
benchmark models we consider here. The units for the 
mass parameters are GeV, and the universal trilinear cou- 
pling A is set to for all models (see |Battaglia et al. 2 003 
for details). 



eration here, which are shown in Figs. |H| and [5] Fig[S] shows 
in particular the differential electron (left) and photon (right) 
yields per neutralino annihilation multiplied by (cru)o, i.e. the 
source function Q(r,E) divided by the number density of neu- 
tralino pairs N P nirs(r) as a function of the particles' kinetic 
energy. As mentioned above, the Bulk and Funnel cases are 
very similar between each other, though in the latter case one 
has a heavier spectrum and a larger value of (av)o- Fig[5]shows 
the same quantity for neutrinos and protons. 

The products of the neutralino annihilation which are more 
relevant to our discussion are secondary electrons and pions. 



Model 


BR(bb) 


BR(t+t-) 


BR(W + W~) 


BR(Z°Z°) 




B' (Bulk) 


74% 


19% 


4% 


0% 


7.8 x 


D' (Coann.) 


21% 


61% 


0% 


0% 


8.9 x 


E' (Focus P.) 


1% 


0% 


90% 


8% 


1.7 x 


K' (Funnel) 


88% 


11% 


0% 


0% 


1.1 X 



Table 2. The branching ratios into various final states for 
the four mSUGRA benchmark models of tab. 2] in the last 



column we also indicate 



in units of cm 3 s . 
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The secondary particles produced by neutralino annihilation 
are subject to various physical mechanisms: i) decay (which 
is especially fast for pions and muons); ii) energy losses which 
can be suffered by stable particles, like electrons and positrons; 
Hi) spatial diffusion of these relativistic particles in the at- 
mosphere of the cluster. Gamma-rays produced by neutral 
pion decay, 7r° — > 77, generate most of the continuum spec- 
trum at energies E ^ 1 GeV and this emission is directly 
radiated since the 7T° — > 77 e.m. decay is very fast. This 
gamma-ray emission is dominant at high energies, ^ 0.3 — 0.5 
of the neutralino mass, but needs to be complemented by 
other two emission mechanisms which produce gamma-rays at 
similar or slightly lower energies: these are the ICS and the 
bremsstrahlung emission by secondary electrons. We will dis- 
cuss the full gamma-ray emission of Coma induced by DM an- 
nihilation in Sect. [I] below. Secondary electrons are produced 
through various prompt generation mechanisms and by the 
decay of charged pions (see, e.g., Colafrancesco & Mele 2001). 
In fact, charged pions decay through ti ± — > /i ± ^ (J (^ (J ), with 
fj, — *• e + + ^ e (^e) and produce e , muons and 

neutrinos. Electrons and positrons are produced abundantly 
by neutralino annihilation (see Fig. |HJ left) and are subject 
to spatial diffusion and energy losses. Both spatial diffusion 
and energy losses contribute to determine the evolution of the 
source spectrum into the equilibrium spectrum of these par- 
ticles, i.e. the quantity which will be used to determine the 
overall multi-wavelength emission induced by DM annihila- 
tion. The secondary electrons eventually produce radiation by 
synchrotron in the magnetized atmosphere of Coma, Inverse 
Compton Scattering of CMB (and other background) photons 
and bremsstrahlung with protons and ions in the atmosphere 
of the Coma cluster (see, e.g.. IColafranccsco & Mele 200l]and 
Colafrancesco 2003 Colafrancesco 2006 for a review). These 
secondary particles also produce heating of the intra-cluster gas 
by Coulomb collisions with the intra-cluster gas particles and 
SZ effect (see, e.g. Colafrancesco 2003 Colafrancesco 2006). 
Other fundamental particles which might have astrophysical 
relevance are also produced in DM annihilation. Protons are 
produced in a smaller quantity with respect to (see Fig. [§] 
right), but do not loose energy appreciably during their lifetime 
while they can diffuse and be stored in the cluster atmosphere. 
These particles can, in principle, produce heating of the intra- 
cluster gas and pp collisions providing, again, a source of sec- 
ondary particles (pions, neutrinos, e , muons, ...) in complete 
analogy with the secondary particle production by neutralino 
annihilation. Neutrinos are also produced in the process of neu- 
tralino annihilation (see Fig.|5] left) and propagate with almost 
no interaction with the matter of the cluster. However, the re- 
sulting flux from Coma is found to be unobservable by current 
experiments. 

To summarize, the secondary products of neutralino an- 
nihilation which have the most relevant astrophysical impact 
onto the multi-frequency spectral energy distribution of DM 
halos are neutral pions and secondary electrons. 



4. Neutralino-induced signals 

A complete description of the emission features induced by DM 
must take, consistently, into account the diffusion and energy- 
loss properties of these secondary particles. These mechanisms 
are taken into account in the following diffusion equation (i.e. 



neglecting convection and re-acceleration effects) 



d dn e 
ffi~dE 



D(E,x)V 



-Q e (E,x) 



dn e 

He 



d_ 

BE 



b(E,x 



dn e 
"dE 



(35) 



where dn e /dE is the equilibrium spectrum, D(E, x) is the dif- 
fusion coefficient, b(E, x) is the energy loss term and Qe(E, x) 
is the source function. The analytical solution of this equa- 
tion for the case of the DM source function is derived in the 
Appendix 1X1 

In the limit in which electrons and positrons lose energy on a 
timescale much shorter than the timescale for spatial diffusion, 
i. e. the regime which applies to the case of galaxy clusters, the 
first term on the r.h.s. of Eq. 1A.11 can be neglected, and the 
expression for equilibrium number density becomes: 



( dE ) nad ^ 



HE) 



dE' Q e (r,E') 



(36) 



(see the Appendix [XJ for a general discussion of the role of 
spatial diffusion and of the regimes in which it is relevant). 

The derivation of the full solution of the diffusion equa- 
tion (Eq. lA.lfl and the effects of diffusion and energy losses 
described in the Appendix [XJ set us in the position to discuss 
the multi-frequency emission produced by the DM (neutralino) 
component of the Coma cluster. We will present the overall 
DM-induced spectral energy distribution (hereafter SED) from 
low to high observing frequencies. 

We describe here our reference setup for the numerical cal- 
culations. Our reference halo setup is the N04 profile and other 
parameters/choice of extrapolation schemes as in Fig. |H| We 
consider the predictions of two particle models, one with a 
branching ratio equal to 1 in bb, i. e. a channel with a soft pro- 
duction spectrum, and the second one with a branching ratio 
equal to 1 into W + W~ , i.e. a channel with hard spectrum. 
Since we have previously shown that diffusion is not relevant 
in a Coma-like cluster of galaxy, we neglect, in our numerical 
calculations, the spatial diffusion for electrons and positrons: 
this is the limit in which the radial dependence and frequency 
dependence can be factorized in the expression for the emis- 
sivity. 

4.1. Radio emission 

At radio frequencies, the DM-induced emission is dominated 
by the synchrotron radiation of the relativistic secondary elec- 
trons and positrons of energy E — 7m e c 2 , living in a mag- 
netic field B(r) and a background plasma with thermal elec- 
tron density n(r), and in the limit of frequency v of the emitted 
photons much larger than the non-relativistic gyro-frequency 
vo — eB/(2nmc) ~ 2.8-B M Hz and the plasma frequency 

v v = 8980 (n(r) /lcm~ 3 ^) 1 2 Hz. Averaging over the directions 
of emission, the spontaneously emitted synchrotron power at 
the frequency v is given by fLongair" 1994 ) : 

f synch (v, E, r) = / d6»^-2yr\/3romc!/osin6'F(a:/sin6') ,(37) 
Jo 2 

where we have introduced the classical electron radius r ( ) = 
e 2 /(me 2 ) = 2.82- 10~ 13 cm, and we have defined the quantities 
x and F as: 



2u 

3^o7 2 



3/2 



(38) 
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and 



F(t) = t dzK 5/3 (z) 



1.25t 1/3 exp( 



-t) [648 + 1 



21 1/12 



(39) 



Folding the synchrotron power with the spectral distribution 
of the equilibrium number density of electrons and positrons, 
we get the local emissivity at the frequency v. 



Jsynch 



dE 



' ' dn p 



dE 



+ 



dn„ 



dE 



- synch 



(40) 



This is the basic quantity we need in order to compare our pre- 
dictions with the available data. In particular, we will compare 
our predictions with measurements of the integrated (over the 
whole Coma radio halo size) flux density spectrum: 



'S'synch(^) — 



,3 jsynch {y, r) 



(41) 



where Dcoma is the luminosity distance of Coma, and with the 
azimuthally averaged surface brightness distribution at a given 
frequency and within a beam of angular size AQ (PSF): 



-^synch 



rffi 



dl 



jsynch {V, I) 



4tt 



(42) 



where the integral is performed along the line of sight (l.o.s.) 
I, within a cone of size Af2 centered in a direction forming an 
angle O with the direction of the Coma center. 

We started from the full dataset on the radio flux density 
spectrum (IThierbach et al. 200311 and minimized the fit with 
respect to the WIMP mass (with the bound M x > 10 GeV for 
the bb case, and mass above threshold for the W + W~ case), 
the strength of the magnetic field (with the bound _B M > 1/iG) 
and the annihilation rate {av)o. The spectrum predicted by 
two models with the lowest values of \r ar e shown in Fig. I1UI 
In both cases the best fit corresponds to the lowest neutralino 
mass allowed, since this is the configuration in which the fall- 
off of the flux density at the highest observed frequency tends 
to be better reproduced. For the same reason, the fit in the 
case of a soft spectrum is favored with respect to the one with 
a hard spectrum (we have checked that in case of t + t~ again 
does not give a bend-over in the spectrum where needed). The 
values of the annihilation rates required by the fit are fairly 
large: {av}o = 4.7 • 10~ 25 cm 3 s _1 for bb case, and about one 
order of magnitude smaller, (av)o — 8.8 ■ 10~ 26 cm 3 s _1 , for 
the W + W~ case, despite the heavier neutralino mass, since the 
best fit values correspond to different values of the magnetic 
field of about 1.2fiG and 8fiG, respectively. 

In Fig. 1111 we compare the radio-halo brightness data 
of |D eiss et al. 19971 with the surface brightness distribution 
I S ynch(r) predicted at v = 1.4 GHz, within a beam equal 
to the detector angular resolution (HPBW of 9.'35), for the 
best fit model with M x = 40 GeV. In the left panel we plot 
the predicted surface brightness considering the case of a uni- 
form magnetic field equal to 1.2/iG, showing explicitly in this 
case that the assumption we made of neglecting spatial diffu- 
sion for electrons and positrons is indeed justified, since the 
results obtained including or neglecting spatial diffusion essen- 
tially coincide. The radial brightness we derive in this case does 
not match the shape of the radio halo indicated by the data. 
However, it is easy to derive a phenomenological setup with a 
magnetic field B(r) varying with radius in which a much bet- 
ter fit can be obtained, while leaving unchanged the total radio 




V [ MHz ] 

Fig. 10. Best fit models for the radio flux density spec- 
trum, in case of a soft spectrum due to a bb annihilation 
final state (solid line, model with M x = 40 GeV) and of 
a hard spectrum due to a W + W~ channel (dashed line, 
model with M x = 81 GeV); values of all parameters set- 
ting the model are given in the text. The datasets is from 
IThierbach et al. 20031 

flux density S a ynch.(v). We show in the right panel of Fig llll 
the predictions for I S y nc h(r) considering a radial dependence 
of the magnetic field of the form: 



B(r) = B 1 



r 

r C 2 



which is observationally driven by the available information on 
the Faraday rotation measures (RM) for Coma (see Fig. 16). 
Such B(r) profile starts at a slightly smaller value in the cen- 
ter of the Coma, rises at a first intermediate scale r c \ and then 
drops rather rapidly at the scale r C 2. The basic information we 
provide here is that a radial dependence of the magnetic field 
like the previous one is required in DM annihilation models to 
reproduce the radio-halo surface brightness distribution. The 
specific case displayed is for best-fit values Bo = 0.55 fj,G , 
P = 2.7, r c i = 3' r C 2 = 17'5, and it provides an excellent fit to 
the surface brightness radial profile (see Fig lllH . In that figure 
we also plot separately the contributions to the surface bright- 
ness due to the smooth DM component (essentially a point-like 
source in case of this rather poor angular resolution) and the 
term due to subhalos (which extends instead to larger radii). 
It is interesting to note that the surface brightness profile can 
only be fitted by considering the extended sub-halo distribu- 
tion which renders the DM profile of Coma more extended than 
the smooth, centrally peaked component. This means that any 
peaked and smooth DM profile is unable to fit this observable 
for Coma. 

A decrease of B(r) at large radii is expected by general con- 
siderations of the structure of radio-halos in clusters and, 
more specifically, for Coma (see Colafrance sco et al. 2 005 1 
and it is also predicted by numerical simulations (see, e.g., 
|Dolag et al. 2002} : thus it seems quite natural and motivated. 
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Fig. 11. Surface brightness distribution at frequency v — 1.4 GHz, within a beam equal to 9'35 (HPBW), for the 
lightest WIMP model displayed in Fig. In the left panel we show the predictions for a model with magnetic field 
that does not change with radius, and in the limit in which spatial diffusion for electrons and positrons has been 
neglected (solid line) or included (dashed line). In the right panel we consider a model with magnetic field taking a 



radial dependence B(r) = B [ 1 + — 



1 + hf- 



with B = 0.55 fiG , /3 = 2.7, r cl = 3' r c2 = 17'.5 in 



order to reproduce the measured surface brightness. In both cases, the contributions from the smooth dark matter 
halo component only and from substructures only are also displayed. The dataset is from lDeiss et al. 19971 



At small radii, the mild central dip of B(r) predicted by the 
previous formula is what is phenomenologically required by the 
specific DM model we worked out in our paper. Finally, we no- 
tice that our specific phenomenological model for the spatial 
distribution of B(r) is able to reproduce the spatial distribu- 
tion of Faraday rotation measures (RMs) observed in Coma 
(see|K im et al. 1990L as shown in Fig |12l It is evident that 
models in which B is either constant or decreases monotoni- 
cally towards large radii seem to be difficult to reconcile with 
the available RM data. The RM data at 9 ^ 20 arcmin seem to 
favour, indeed, a model for B(r) with a slight rise at interme- 
diate angular scales followed by a decrease at large scales, like 
the one we adopt here to fit the radio-halo surface brightness 
of Coma. In this respect, it seems that our choice for B(r) is, 
at least, an observationally driven result. 

The synchrotron signal produced by the annihilation of 
DM depends, given the fundamental physics and astrophysics 
framework, on two relevant quantities: the annihilation rate 
and the magnetic field. Thus, it is interesting to find the best- 
fitting region of the {av)o — B plane which is consistent with 
the available dataset for Coma. 

Since the data on the radio flux density spectrum 
( Thicrbac h et al. 200311 is a compilation of measurements per- 
formed with different instruments (possibly with different sys- 
tematics), it is difficult to decide a cut on the Xr value which 
defines an acceptable fit. In Fig. 1131 we plot sample isolevel 
curves for Xr> spotting the shape of the minima of Xr i i n the 
plane B — {uv)o, for the two sample annihilation channels and 
a few sample values of the WIMP mass (note that values label- 
ing isolevel curves are sensibly different in the two panels). In 
Fig. I14l we show the analogous Xr isolevels in the WIMP mass 
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Fig. 12. The observed absolute RM of background sources in 
the Coma field are shown as a function of projected radius 6 
in arcmin. The blu curve is the prediction of a model with 
B — const. The red dashed curve is the prediction of a model 
with B(r) decreasing monotonically towards large radii and the 
solid red curve is the prediction of our model that fits the Coma 
radio-halo surface brightness. Data on positive (filled dots) and 
negative (empty dots) RMs are from lKim et al. 19901 



- annihilation rate plane, and taking at each point the min- 
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imum Xr while varying the magnetic field strength between 
1/iG and 20fj,G: the curves converge to a maximal value en- 
forced by the lower limit of lfJ,G, and the upper value does not 
enter in defining isolevel curve shapes. 

In order to assess whether the outlined radio-data preferred 
regions are or are not compatible with supersymmetric DM 
models, we proceed to a random scan of the SUSY parame- 
ter space, in the bottom-up approach which we outline below. 
We relax all universality assumptions, and fully scan the low- 
energy scale MSSM, imposing phenomenological as well as cos- 
mological constraints on the randomly generated models 4 . We 
take values of tan /3, the ratio between the vacuum expecta- 
tion values of the two Higgs doublets, between 1 and 60. The 
parameters entering the neutralino mass matrix are generated 
in the range 1 GeV < mi, m2, fJ, < 1 TeV, and we define 
jtilsp = min(mi, rri2, fJ,). To avoid flavor changing effects in 
the first two lightest quark generations, we assume that the 
soft-breaking masses in the first two generations squark sec- 
tor are degenerate, i.e. we assume mli' ~ ~ = m~'~ ~. The 

Q, U, D Q, U, D 



WMAP range, and with red circles those producing a relic 
abundance below the WMAP range (whose relic abundance 
can however be cosmologically enhanced, in the context of 
quintessential or Brans-Dicke cosmologies, or which can be 
non-thermally produced, as to make up all of the observed 
CDM (see IMurakami fc Wells 2001 and other refe.V The low 
\ 2 ranges of {av)o — B and (av)o — M x values indicated in 
Figs. [TJl and 1141 are therefore shown to be actually populated 
by a number of viable SUSY models. 



scalar masses are scanned over the range 



m-LSP < m 



(1 < 3) mL M) , 
u 



,(1,3) 



(1,2,3) (1,2,3) 



4.2. From the UV to the gamma-ray band 

Inverse Compton (IC) scatterings of relativistic electrons and 
positrons on target cosmic microwave background (CMB) pho- 
tons give rise to a spectrum of photons stretching from below 
the extreme ultra-violet up to the soft gamma-ray band, peak- 
ing in the soft X-ray energy band. Let E = r ym e c 2 be the 
energy of electrons and positrons, e that of the target pho- 
tons and E 1 the energy of the scattered photon. The Inverse 
Compton power is obtained by folding the differential number 

density of target photons with the IC scattering cross section: 
m A < 2.5TeV(43) 
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Pic (E-y,E) = cE 1 I den(e)a(E^,e,E) 



(49) 



where n(e) is the black body spectrum of the 2.73-R' CMB pho- 
tons, while a(E-,, e, E) is given by the Klein-Nishina formula: 



<7(£ 7 ,e,£) 
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4e7 2 



G(q,T e 



Finally, we take the gluino mass in the range 200 GeV < 
m~ < 3 TeV. The mass ranges for squarks and gluino have 
been chosen following qualitative criteria (IBaer et al. 20031 
|Battaglia et al. 2003| l, so that all viable models generated 
should be "visible" at the LHC. 

We exclude models giving a relic abundance of neutralinos 
exceeding Q x h 2 > 0.13. Further, we impose the various collid- 
ers mass limits on charginos, gluinos, squarks and sleptons, as 
well as on the Higgs masses 5 . Moreover, we also require the 
BR(& — > s'y) and all electroweak precision observables to be 
consistent with the theoretical and experimental state-of-the- 
art JEidelman et al. 20041 . 

We classify the models according to the branching ratios of 
the neutralino pair-annihilations final states, according to the 
following criteria: we consider a model having a hard spectrum 
if 



where or is the Thomson cross section and 

(r e g) 2 (i-g) 



G(q,T 
with 



2qlnq + (l + 2q)(l-q) + 



2(i + r e(? ) 



T e = 4e 7 /(m e c 2 ) q = £ 7 / [r e ('ym e c 2 - E-,)] 



(50) 



(51) 



(52) 



Folding the IC power with the spectral distribution of the equi- 
librium number density of electrons and positrons, we get the 
local emissivity of IC photons of energy E*y : 
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(53) 



which we use to estimate the integrated flux density spectrum: 



BR(VU+VU~) + BR{ZZ) > 0.; 



(47) Si C (£ 7 ) = / d 3 r 



. he (Ej 



(54) 



a soft spectrum is instead attributed to models satisfying the 
condition 



6 

E 

i=l 



[BR(M0 + BR(g i g iS )] + BR(j 9 ) > 0.8. 



(48) 



We show, in Fig |15l a scatter plots of the viable SUSY con- 
figurations, indicating with filled green circles those ther- 
mally producing a neutralino relic abundance within the 2-cr 



4 We scan all the SUSY parameters linearly over the indi- 
cated range. 

5 Since we do not impose any gaugino unification relation, 
we do not impose any constraint from collider searches on the 
neutralino sector. 



In Eq. 1491 and Eq. 153H the limits of integration over e and E-, 
are set from the kinematics of the IC scattering which restricts 
q in the range l/(47 2 ) < q < 1. 

The last relevant contribution to the photon emission of 
Coma due to relativistic electrons and positrons is the process 
of non-thermal bremsstrahlung, i. e. the emission of gamma-ray 
photons in the deflection of the charged particles by the electro- 
static potential of intra-cluster gas. Labeling with E — r ym e c 2 
the energy of electrons and positrons, and with _E 7 the energy 
of the emitted photons, the local non-thermal bremsstrahlung 
power is given by: 



Pb{E-,,E, 
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Fig. 13. Isolevel curves for reduced \ 2 from the fit of the Coma radio flux density spectrum, for a given mass (50 
- 100 -150 GeV) and annihilation channel. The halo profile is the best fit N04 profile: M v „ = 0.9 10 15 M Q hr 1 and 
c v ir — 10, with subhalo setup as in Fig. El 




Fig. 14. Isolevel curves for minimum reduced \ 2 from the fit of the Coma radio flux density spectrum, obtained by 
varying the magnetic field within 1//G < < 20/^G, for a given annihilation channel. The halo profile is the best fit 
N04 profile: M vir = 0.9 1O 15 M h^ 1 and c vir = 10, with subhalo setup as in Fig.© 



with the sum including all species j in the intra-cluster 
medium, each with number density rij(r) and relative produc- 
tion cross section: 



<Ti(Bh,E) 



8nE~, 



(l + (l-E-,/Ef 



^ - |(1 



Ej/E) i 



where a is the fine structure constant, <f>i and <j>2 two energy de- 
pendent scattering functions which depend on the species j (see 
|Longair 1994[ for details). The emissivity Jb {E 1 ,r) is obtained 
by folding the power over the equilibrium electron/positron 
number density, i.e. the analogous of Eq. 15311 . while the in- 
tegrated flux density 5b (-E7) is obtained by summing over all 
relevant sources as in Eq. 1541 . We apply this scheme to Coma 
implementing the gas density profile in Eq. by including 



atomic and molecular hydrogen and correcting for the helium 
component. 

As we have already mentioned, a hard gamma-ray com- 
ponent arises also from prompt emission in WIMP pair anni- 
(Sfijations, either in loop suppressed two-body final states giv- 
ing monochromatic photons, or through the production and 
prompt decay of neutral pions giving gamma-rays with contin- 
uous spectrum. Since photons propagate on straight lines (or 
actually geodesies), the gamma-ray flux due to prompt emis- 
sion is just obtained by summing over sources along the line of 
sight; we will consider terms integrated over volume 



F 7 (£ 7 ) 
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Fig. 15. A scatter plot of SUSY models, consistent with all available phenomenological constraints, giving a relic 
abundance in the 2-er WMAP range (green filled circles) or below it (low, relic density models, red circles), for soft 
(left panel) and hard (right panel) annihilation spectra. 
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Fig. 16. Multi- wavelength spectrum of the two best fit models for the radio flux shown in Fig. (see text for details). 
The halo profile is the best fit N04 profile: M V i r = 0.9 10 15 M© and c V i r = 10, with subhalo setup as in Fig. El 



4.3. The multi-frequency SED of Coma 

We show in Fig. 1161 the multi-frequency SED produced by 
WIMP annihilation in the two models used to fit the radio 
halo spectrum of Coma, as shown in Fig. I1UI 

The model with M x = 40 GeV provides the better fit 
to the radio halo data because the relative equilibrium elec- 
tron spectrum is steeper and shows also the high-f bending 
which fits the most recent data (Thicrba ch et al. 2003t . The 
IC and bremsstrahlung branches of the SED are closely related 
to the synchrotron branch (since they depend on the same par- 
ticle population) and their intensity ratio depends basically on 
the value of the adopted magnetic field. The relatively high 
value B = 1.2 /j,G indicated by the best fit to the radio data 
implies a rather low intensity of the IC and bremsstrahlung 
emission, well below the EUV and hard X-ray data for Coma. 



Nonetheless, the gamma-ray emission due to 7T 77 decay 
predicted by this model could be detectable with the GLAST- 
LAT detector, even though it is well below the EGRET upper 
limit. 

The detectability of the multi-frequency SED worsens in 
the model with M x = 81 GeV, where the flatness of the equi- 
librium electron spectrum cannot provide an acceptable fit 
to the radio data. Moreover, the adopted value of the mag- 
netic field B — 8 fiG implies a very low intensity of the IC, 
bremsstrahlung and ty° — > 77 emission, which should be not 
detectable by the next generation HXR and gamma-ray exper- 
iments. 

The energetic electrons and positrons produced by WIMP 
annihilation have other interesting astrophysical effects among 
which we will discuss specifically in the following the Sunyaev- 
Zel'dovich (hereafter SZ) effect produced by DM annihilation 
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and the heating of the intracluster gas produced by Coulomb 
collisions. 



4.4. SZ effect 

The energetic electrons and positrons produced by WIMP an- 
nihilation interact with the CMB photons and up-scatter them 
to higher frequencies producing a peculiar SZ effect (as origi- 
nally realized by Colafrancesco 2004 1 with specific spectral and 
spatial features. 

The generalized expression for the SZ effect which is valid in 
the Thomson limit for a generic electron population in the rela- 
tivistic limit and includes also the effects of multiple scatterings 
and the combination with other electron population in the clus- 
ter atmospheres has been derived by Colafranccsco ~ al. 20031 
This approach is the one that should be properly used to calcu- 
late the specific SZdm effect induced by the secondary electrons 
produced by WIMP annihilation. Here we do not repeat the 
description of the analytical technique and we refer to the gen- 
eral analysis described in ICo lafranccsco et al. 2003 According 
to these results, the DM induced spectral distortion is 



A7 D m(s) = 2 



(fc B r ) 

(he) 2 



-Vom g(x) 



(58) 



where To is the CMB temperature and the Comptonization 
parameter i/dm is given by 



J/DM 



(It 



(59) 



in terms of the pressure Pdm contributed by the secondary 
electrons produced by neutralino annihilation. The quantity 
2/dm oc (av)on\ and scales as oc (av) M~ 2 , providing an in- 
creasing pressure Pdm and optical depth tdm = ct / d£n e 
for decreasing values of the neutralino mass M x . The function 
g(x), with x = hv/kBTo, can be written as 



g(x) = 



(k B T c ) 



io(xe a )P(s)ds — io(x) 



in terms of the photon redistribution function P(s) and of 
x 3 /(e x — 1), where we defined the 



io(x) = 2(k B T ) 3 /(hc) 2 
quantity 



<fc B T c > = ^ 
T 



Pdl = 



f n c d£ 



dpf e (p)-pv(p)m c c (61) 



(see IColafrancc sco et al. 20031 IColafrancesco "20041 . which is 
the analogous of the average temperature for a thermal pop- 
ulation (for a thermal electron distribution (ksTc) = k B T c 
obtains, in fact). The photon redistribution function P(s) = 
J dpf c (p)P B (s;p) with s = ln(z//V), in terms of the CMB pho- 
ton frequency increase factor v' /v = I7 2 — |, depends on 
the electron momentum (p) distribution, f c (p), produced by 
WIMP annihilation. 

We show in Fig |17l the frequency dependence of the CMB 
temperature change, 



AT 



(e x 



l) 2 Al 



(62) 



as produced by the DM-induced SZ effect in the two best fit 
WIMP models here considered, compared to the temperature 
change due to the thermal SZ effect produced by the intraclus- 
ter gas. The most recent analysis of the thermal SZ effect in 




100 200 300 400 500 600 700 800 900 1000 
frequency (GHz) 



Fig. 17. The SZ effect produced by the bb model with 
M x = 40 GcV (black solid curve) and by the W + W~ 
model with M x — 81 GeV (black dashed curve) in Coma 
are shown in comparison with the thermal SZ effect of 
Coma (blue curve). The red curves represent the overall 
SZ effect. Notice that the DM-induced SZ effect has a 
very different spectral behavior with respect to the ther- 
mal SZ effect. SZ data are from OVRO (magenta), WMAP 
(cyan) and MITO (blue). The sensitivity of PLANCK (18 
months, 1 a) is shown for the LFI detector at 31.5 and 
53 GHz channels (cyan shaded regions) and for the HFI 
detector 143 and 217 GHZ channels (green and yellow 
shaded areas, respectively). 



Coma (DcPct ris et al. 200311 provides an estimate of the opti- 
(60) cal depth of the thermal intracluster gas rth = 4.9- 10 -3 which 



a detectable SZdm effect which has a quite different spectral 
shape with respect to the thermal SZ effect: it yields a tem- 
perature decrement at all the microwave frequencies, ^$ 600 
GHz, where the thermal SZ effect is observed and produces a 
temperature increase only at very high frequencies > 600 GHz. 
This behavior is produced by the large frequency shift of CMB 
photons induced by the relativistic secondary electrons gener- 
ated by the WIMP annihilation. As a consequence, the zero 
of the SZdm effect is effectively removed from the microwave 
range and shifted to a quite high frequency ~ 600 GHz with 
respect to the zero of the thermal SZ effect, a result which al- 
lows one, in principle, to estimate directly the pressure of the 
electron populations and hence to derive constraints on the 
WIMP model (see IColafrancesco 200-fl . 

The presence of a substantial SZdm effect is likely to dom- 
inate the overall SZ signal at frequencies x <^ 3.8 — 4.5 pro- 
viding a negative total SZ effect. It is, however, necessary to 
stress that in such frequency range there are other possible 
contributions to the SZ effect, like the kinematic effect and 
the non-thermal effect which could provide additional biases 
(see, e.g., Colafrancesco et al. 20031. Nonetheless, the peculiar 
spectral shape of the SZum effect is quite different from that of 
the kinematic SZ effect and of the thermal SZ effect and this 
result allows us to disentangle it from the overall SZ signal. 
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Fig. 18. The ratio between the DM-induced and the ther- 
mal SZ effect in Coma is shown for the model with 
M x = 40 GeV (red curve) and for the model with M x = 81 
GeV (blue curve). The model with M x = 40 GeV produces 
an SZ effect which could be detected with the next coming 
microwave experiments. 

An appropriate multi-frequency analysis of the overall SZ ef- 
fect based on observations performed on a wide spectral range 
(from the radio to the sub-mm region) is required to separate 
the various SZ contributions and to provide an estimate of the 
DM induced SZ effect. In fact, simultaneous SZ observations at 
low frequencies ~ 30 GHz (where there is the largest tempera- 
ture decrement due to SZdm), at ~ 150 GHz (where the SZdm 
deepens the minimum in A I /I with respect to the dominant 
thermal SZ effect), at ~ 220 GHz (where the SZdm dominates 
the overall SZ effect and produces a negative signal instead of 
the expected « null signal) and at 250 GHz (where the still 
negative SZdm decreases the overall SZ effect with respect to 
the dominant thermal SZ effect) coupled with X-ray observa- 
tions which determine the gas distribution within the cluster 
(and hence the associated dominant thermal SZ effect) can sep- 
arate the SZdm from the overall SZ signal, and consequently, 
set constraints on the WIMP model. 

The WIMP model with M x = 40 GeV produces a temper- 
ature decrement which is of the order of ~ 40 to 15 piK for 
SZ observations in the frequency range ~ 30 to 150 GHz (see 
Fig |17t , These signals are still within the actual uncertainties 
of the available SZ data for Coma and are below the current 
SZ sensitivity of WMAP (see, e.g,. IBennet et al. 20031 and the 
results of the analysis of the WMAP SZ signals from a sample 
of nearby clusters performed bv lLieu et al. 2005 1 . Nonetheless, 
such SZ signals could be detectable with higher sensitivity ex- 
periments. The high sensitivity planned for the future SZ ex- 
periments can provide much stringent limits to the additional 
SZ effect induced by DM annihilation. In this context, the next 
coming sensitive bolometer arrays (e.g., APEX), interferomet- 
ric arrays (e.g., ALMA) and the PLANCK-HFI experiment, or 
the planned OLIMPO balloon-borne experiment, have enough 
sensitivity to probe the contributions of various SZ effects in 
the frequency range ~ 30 — 250 GHz, provided that accu- 
rate cross-calibration at different frequencies can be obtained. 



The illustrative comparison (see Fig |17t between the model 
predictions and the sensitivity of the PLANCK LFI and HFI 
detectors at the optimal observing frequencies (y = 31.5 and 
53 GHz for the LFI detector and v = 143 and 217 GHz for 
the HFI detector) show that the study of the SZ effect pro- 
duced by DM annihilation is actually feasible with the next 
generation SZ experiments. We show in Fig |18l the expected 
ratio between the DM-induced SZ effect and the thermal SZ 
effect for the two WIMP models here considered. It is evident 
that while the model with M x = 40 GeV provides a detectable 
signal which is a sensitive fraction of the thermal SZ effect 
at v < 250 GHz, the SZ signal provided by the model with 
M x = 81 GeV is by far too small to be detectable at any fre- 
quency. 

The spectral properties shown by the SZdm for neutralinos 
depends on the specific neutralino model as we have shown in 
Fig |17l in fact, the SZ effect is visible for a neutralino with 
M x = 40 GeV and not visible for a neutralino with M x =81 
GeV. Thus the detailed features of the SZ effect from DM anni- 
hilation depends strongly on the mass and composition of the 
DM particle, and - in turn - on the equilibrium spectrum of the 
secondary electrons. Each specific DM model predicts its own 
spectrum of secondary electrons and this influences the relative 
SZ effect. Models of DM which provide similar electron spectra 
will provide similar SZ effects. 

4.5. Heating of the intracluster gas 

Low energy secondary electrons produced by WIMP annihi- 
lation might heat the intracluster gas by Coulomb collisions 
since the Coulomb loss term dominates the energy losses at 
E ^ 200 MeV (see Fig. IA.3I . The specific heating rate is given 
by 

dE f dn e (dE\ 

where ^nf- is the equilibrium electron spectrum derived 
in Sect. ^ and the Coulomb loss rate is (dE/dt)coui = 
b° Coul n (1 + log(7/n)/75) where n is the mean number den- 
sity of thermal electrons in cm' " 3 (see Eq. EU the average 
over space gives about n ~ 1.3 10 -3 ), 7 = E/m e and 
b° Coul ~ 6.13- 10" 16 GeVs" 1 . Fig[T9lshows the specific heating 
rate of Coma as produced in the two WIMP models explored 
here. The non-singular N04 halo model adopted in our analy- 
sis does not provide a high specific heating rate at the cluster 
center, and thus one might expect an overall heating rate for 
Coma which is of order of ~ 10 38 erg/s (~ 10 36 erg/s) for 
the WIMP model with M x = 40 GeV (M x = 81 GeV). We 
also notice that the region that mostly contributes to the over- 
all heating of Coma is not located at the center of the cluster. 
This is again a consequence of the non-singular N04 DM profile 
which has been adopted. The diffusion of electrons in the in- 
nermost regions of Coma acts in the same direction and moves 
the maximum of the curves shown in the right panel of Fig |19l 
towards the outskirts of Coma, even in the case of a halo den- 
sity profile which is steeper than the adopted one. 
This implies, in conclusion, that WIMP annihilation cannot 
provide most of the heating of Coma, even in its innermost 
regions. Such a conclusion seems quite general and implies 
that non-singular DM halo models are not able to provide 
large quantities of heating at the center of galaxy clusters so 
to quench efficiently the cooling of the intracluster gas (with 
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Fig. 19. Left: the specific heating rate is plotted against 
the radial distance from the center of Coma. Right: the 
specific heating rate multiplied by the volume element 
is plotted against the radial distance from the center of 
Coma. 

powers of ~ iq 43-45 erg/s). Only very steep halo profiles (even 
steeper than the Moore profile) and with the possible adia- 
batic growth of a central matter concentration (e.g., a central 
BH) could provide sufficient power to quench locally (i.e. in 
the innermost regions) the intracluster gas cooling (see, e.g., 
Totani 2004). However, we stress that the spatial diffusion of 
the secondary electrons in the innermost regions of galaxy clus- 
ters should flatten the specific heating rate in the vicinity of 
the DM spike and thus decrease substantially the heating effi- 
ciency by Coulomb collisions. In conclusion, we believe that the 
possibility to solve the cooling flow problem of galaxy clusters 
by WIMP annihilation is still an open problem. 

5. Discussion 

WIMP annihilation in galaxy cluster is an efficient mecha- 
nism to produce relativistic electrons and high-energy particles 
which are able, in turn, to produce a wide SED extended over 
more than 18 orders of magnitude in frequency, from radio to 
gamma-rays. We discuss here the predictions of two specific 
models which embrace a vast range of possibilities. 

The bb model with M x = 40 GeV and annihilation cross 
section {av)o = 4.7 ■ 10 _25 cm 3 s _1 provides a reasonable fit to 
the radio data (both the total spectrum and the surface bright- 
ness radial distribution) with a magnetic field whose mean 
value is B m 1.2/^G. We remind here that the quite high value 
of {crv)o is well inside the range of neutralino masses and an- 
nihilation cross-sections provided by the most general super- 
symmetric DM setup (see our discussion in Sect. 3. 2). Table [3] 
provides an illustrative scheme of the radiation mechanisms, of 
the particle energies and of the fluxes predicted by this best-fit 
WIMP model for a wide range of the physical conditions in the 
cluster atmosphere. 



For the best-fit values of M x — 40 GeV and (av)o = 
4.7- 10~ 25 cm 3 s _1 this model yields EUV and HXR fluxes which 
are more than one order of magnitude fainter than the Coma 
data. The gamma-ray flux produced by this model is domi- 
nated by the continuum ir° — > 77 component and it is a fac- 
tor ~ 5 lower than the EGRET upper limit of Coma at its 
peak frequency (see Fig. 1161 left panel). Such gamma-ray flux 
could be, nonetheless, detectable by the GLAST-LAT detec- 
tor (we will discuss more specifically the detectability of the 
gamma-ray WIMP annihilation signals from galaxy clusters 
in a dedicated, forthcoming paper (Colafrancesco, Profumo & 
Ullio 2006b). The rather low neutralino mass M x = 40 GeV of 
this model makes it rather difficult to be testable by Cherenkov 
gamma-ray detectors operating at higher threshold energies. 

Increasing the neutralino mass does not provide a good fit 
of the radio-halo spectrum (see Fig. 1161 right panel) and yields, 
in addition, extremely faint EUV, HXR and gamma-ray fluxes, 
which turn out to be undetectable even by GLAST and/or by 
the next coming high-energy experiments. 

It is possible to recover the EUV and HXR data on Coma 
with a IC flux by secondary electrons by increasing the an- 
nihilation cross-sections by a factor ~ 10 2 (i.e., up to values 
(<tv) » 7 • 10~ 23 cm 3 s _1 ) in the best-fit bb soft WIMP model 
(at fixed M x = 40 GeV). However, in such a case both the 
radio-halo flux and the hard gamma-ray flux at ~ 1 GeV as 
produced by n° decay should increase by the same factor lead- 
ing to a problematic picture: in fact, while the radio-halo data 
would imply lower values of the magnetic field B ~ 0.1 /iG 
which might still be allowed by the data, the ir° — > 77 gamma- 
ray flux at E > 100 MeV should exceed the EGRET limit on 
Coma. This option is therefore excluded by the available data. 
Alternatively, it would be possible to fit the EUV and HXR 
spectra of Coma with the adopted value of (crv)o ~ 7 ■ 
10" 23 cm 3 s~ 1 for the bb model with M x = 40 GeV, in the case 
we sensibly lower the mean magnetic field. Values of the aver- 
age magnetic field ^ 0.2/iG are required to fit the HXR flux of 
Coma under the constraint to fit at the same time the radio- 
halo spectrum (see Fig |20l left panel), consistently with the 
general description of the ratio between the synchrotron and IC 
emission powers in Coma (see, e.a..lColafrance sco et al.' 2005 
Rci mer et al. 20 04). Magnetic fields as low as ~ 0.15pG can fit 
both the HXR and the EUV fluxes of Coma. However, also in 
this case the 7r° — » 77 gamma-ray flux predicted by the same 
model at E > 100 MeV exceeds the EGRET limit on Coma, 
rendering untenable this alternative. Actually, the EGRET up- 
per limit on Coma set a strong constraint on the combination 
of values B and (crv)o (see Fig |20l right panel) so that mag- 
netic field larger than ^ 0.3 /iG are required for the parameter 
setup of the bb model with M x = 40 GeV. Fig|2Ul shows the 
upper limits on the value of (crv) as a function of the assumed 
value of the mean magnetic field of Coma. According to these 
results, it is impossible to fit all the available data on Coma 
for a consistent choice of the DM model and of the cluster 
magnetic field. The EUV and HXR data in particular require 
extreme conditions, i.e. low values of the magnetic field and/or 
high values of the annihilation cross section, which violate the 
EGRET gamma-ray limit. Thus, realistic DM models that are 
consistent with the radio and gamma-ray constraints predict 
IC emission which falls short of fitting the EUV and HXR data 
of Coma. 

An appealing property of the WIMP model worked out 
here is that it can reproduce both the spatial distribution 
of the radio-halo surface brightness of Coma and, in princi- 
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Table 3. Predicted flux in various frequency ranges for a neutralino bb model with M x = 40 GeV and (av) 
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pie, also the spatial profile of the EUV emission (see, e.g., 
|Bowyer et al. 2004^ which seems more concentrated than the 
radio-halo surface brightness. As for the radio-halo surface 
brightness profile, it seems necessary for this WIMP model 
- due to the shape of the DM halo profile - to invoke a radial 
distribution of the magnetic field with a mild decrease towards 
the Coma center to counterbalance the centrally peaked DM 
profile, and with an exponential cutoff at large radii to counter- 
balance the effect of the subhalo distribution. We notice here 
that such a specific B(r) spatial distribution is - interestingly 
enough - able to reproduce the radial distribution of the RMs 
in Coma (see Sect. 5.1). While the (Synchrotron) radio surface 
brightness depends strongly on the magnetic field radial pro- 
file, the (ICS) EUV radial profile only depends on the DM 
halo profile and on the secondary electron properties and is, 
hence, more concentrated (see Fig. 1101 for an example). Thus, 
the radial distribution of the EUV emission could be reason- 
ably reproduced by the WIMP model which best fits the radio 
data but with a very low value of the magnetic field of order of 
;$ 0.15 fiG. We already noticed, however, that models with val- 
ues of the average magnetic field in Coma which are ^ 0.3/iG 
produce a gamma-ray flux which exceeds the EGRET upper 
limit of Coma (see Fig. I2()|l . rendering these models untenable. 

We summarize all the constraints on the neutralino models 
set by the magnetic field and by the annihilation cross-section 
in Figs. 1201 and 1211 The available data set constraints on the 
WIMP annihilation rate. Figs. 22 and 23 show the upper lim- 
its on {av)o as a function of the assumed value for the mean 
magnetic field in Coma. The EGRET limit proves to be the 
more constraining at the moment with respect to the HXR 
and EUV data. These limits are able to test directly the anni- 
hilation rate since they are independent of the magnetic field 
value. Nonetheless, the combination of the gamma-ray and/or 
HXR constraints with the radio constraints will be able to 
determine the full setup of the relevant quantities whose com- 
bination is able to fit the overall Coma SED. In this context, it 
is clear that the possible GLAST observations of Coma, com- 
bined with the radio data, will increase by far the constraints 
in the (av)o — B plane. 

The results of our analysis also depend on the assumed DM 
halo density profile. Fig. l22l shows the scaling of fluxes with the 
assumptions on the halo model for Coma. We compare here, 
for the sake of illustration, the scalings of the N04 and of the 
Burkert model. Switching to one of the halo models displayed 
here is equivalent to shifting all values of {av)o plotted in the 
figures to {av)o divided by the scaling value shown here. This 
analysis allows us to compare correctly the results of the multi- 
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Fig. 23. Multi-wavelength spectra for the four benchmark 
models described in the text. The prediction is shown for 
the best fit N04 profile, and our reference choice for sub- 
halo parameters, and for a mean magnetic field equal to 
2^G. 



frequency analysis we have presented in this paper in terms of 
substructure enhancement and halo density profile. 

Table [I] also shows the typical values for the annihilation 
cross-section (av)o and the relative signals expected at dif- 
ferent frequencies for the Coma best-fit model (bb neutralino 
model with M x = 40 GeV and B = 1 /j,G) that we explored 
in this paper. Large neutralino pair annihilation cross section 
will, in general, produce sizable signals also for other indi- 
rect detection techniques, including antimatter searches and 
gamma rays from the center of the Milky Way. Antimatter and 
gamma-ray fluxes also largely depend on the Milky Way dark- 
matter halo and on the specific neutralino model (e.g. through 
the antimatter yield per neutralino annihilation). Existing 
analysis JBaer and Pro fumo 2005 Profu mo and Ullio 2 004 
Bac r et al. 2005al IProfum~ 2005 1 make it possible to draw 
some qualitative estimates of the cross-section, {crv)o, required 
to produce sensible signals at future DM search experiments. 
We provide in Table [I] order-of-magnitude estimates for the 
value of {av)o expected to give observable signals in space- 
based antimatter (AMS-02, Pamela, GAPS) and gamma-ray 
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Fig. 20. Scaling of the multi-wavelength spectrum and of relative bounds on the particle physics model with the 
assumed value for the mean magnetic field in Coma. Left panel: we have chosen a few sample values for the magnetic 
field and varied freely pair annihilation cross section and WIMP mass to minimize the x 2 f° r the fit of radio data (a 
bb final state is assumed); the decrease in the magnetic field must be compensated by going to larger M x and (crv)o, 
with a net increase in (av)o/M^., as it can be seen from the increase in the 7r° component. The increase in the IC 
component is, at large values for the magnetic field, significantly more rapid, since for large values of the magnetic field 
synchrotron losses are the main energy loss mechanism for electrons and positrons and tend to decrease the number 
density for the equilibrium population. Right panel: upper limit on (crv)o as a function of the assumed value for the 
mean magnetic field in Coma; at each wavelength the limit is derived assuming that the predicted flux should be lower 
than the upper limit from each individual data point (slight overestimate of the limit from radio data, but we do not 
need to decide the cut on the reduced % 2 marking the overshooting of the radio flux) . Two sample values of M x are 
assumed. The lines marked GLAST refer to the GLAST projected sensitivity assuming no other 7-ray component is 
present. In both panels the halo profile is the best fit N04 profile: M„ r = 0.9 10 15 M Q and c V i r = 10, with subhalo 
setup as in Fig. El 




Fig. 21. The analogous of Fig. |201 but now for the W + W~ final state and fixing the WIMP mass to 81 GeV. In the 
left panel, for each value of the magnetic field, the value of (crw)o is obtained by normalizing the radio flux at the value 
of the flux at the highest frequency point in the available dataset. Note again that on fair fit of the full radio dataset 
can be derived for this hard channel. 



search experiments (GLAST) , for two extreme choices of the galactic dark-matter halo: a cuspy profile (such as the N04 pro- 
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Fig. 22. Scaling of fluxes with the assumptions on the halo model for Coma. In the plane c V i r - A 2 we plot isolevel 
curves for fluxes normalized to the corresponding values within the setup as for the N04 profile in Fig. El marked 
with a dot in the left panel, i.e. the halo model we have assumed so far as reference model; the model marked with 
square in the right panel corresponds to the Burkert profile selected in Fig. For all examples displayed we have 
fixed M V i r = 0.9 10 15 M Q h~ x and a 50% mass in substructures. The left panel refers to N04 profiles, the right panel 
to Burkert profiles; reference values for c V i r and A 2 , as obtained by fitting the corresponding halo profiles and from 
our discussion on the substructure role, are marked on the axis. Switching to one of the halo models displayed here is 
equivalent to shifting all values of (av) plotted in figures to (av) divided by the scaling value shown here. 



file) and a cored profile (such as the Burkert profile). Different 
search techniques, such as direct detection, or neutrino flux de- 
tection from the core of the Sun induced by the annihilation 
of captured neutralinos, critically depend upon the scattering 
cross section of neutralinos off nucleons, and the resulting de- 
tection rates are therefore unrelated, in general, to the pair 
annihilation process discussed here. 

We show in Fig. 1231 the overall SED of Coma as expected 
from the four benchmark models described in the Sect. 3. 2. The 
predictions are shown for the best fit N04 profile and for our 
reference choice for subhalo parameters, and for a mean mag- 
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Table 4. Order-of-magnitude estimates of the values of 
the cross-section, (<jv)q, in units of cm 3 s , needed to 
reproduce the detected non-thermal emission features in 
Coma and to produce sizable signals for future indirect 
dark matter search experiments in the Milky Way. 



netic field of 2fiG. None of these benchmark model may provide 
a reasonable fit to the radio data. Notice, in addition, the quite 
dim multi-frequency SED predicted for Coma in these bench- 
mark models. The largest fluxes are, not surprisingly, obtained 
for the model lying in the focus point region (E'~). In that 
case, neutralinos mostly annihilate into gauge bosons, as can 
be inferred from the spectral shape, which closely resembles 
that in Fig. 1211 We therefore conclude that the expectation 
of astrophysical signatures from neutralino DM annihilations 
in the Coma cluster (with natural assumptions on the dark 
halo profile, substructures and magnetic field of Coma) is not 
promising in the context of the commonly discussed minimal 
supergravity scenario. 

It would be interesting to compare the predictions of the 
WIMP annihilation for Coma with the implication of the pres- 
ence of another population of cosmic rays of different origin 
like that, often invoked, produced by acceleration processes in 
the atmosphere of Coma. Acceleration scenarios usually pro- 
duce power-law spectra for the electrons which are primarily 
accelerated by shocks or turbulence and are, hence remarkably 
different from the source spectra produced by neutralino ac- 
celeration. Specifically, acceleration models do not exhibit a 
cut-off at the neutralino mass and do not produce the pecu- 
liar peaked it gamma-ray emission which remains a distinctive 
feature of neutralino DM models. A continuum ir° gamma-ray 
emission can be produced in secondary models where the elec- 
trons are produced by proton-proton collision in the cluster 
atmosphere. But even in this case the gamma-ray spectrum is 
likely to be resembled by a power-law shape which keeps mem- 
ory of the original acceleration events for the hadrons. Thus, 
it will be possible to separate DM annihilation models from 
acceleration models based on multi-frequency observations of 
the hadronic and leptonic components of the cluster SED. 
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Finally, It should be noticed that the same problem with 
the consistent fitting of both the synchrotron and IC compo- 
nents to the radio and EUV/HXR data of Coma still remains 
in both DM and acceleration models, pointing to the fact that 
these spectral features, if real, have probably different physical 
origin. 

6. Summary and conclusions 

WIMP annihilations in galaxy cluster inevitably produce high- 
energy secondary particles which are able, in turn, to produce 
a wide SED extended over more than 18 orders of magnitude 
in frequency, from radio to gamma-rays. 

A consistent analysis of the DM distribution and of its 
annihilation in the Coma cluster shows that WIMP annihi- 
lation is able to reproduce both the spectral and the spatial 
features of the Coma radio halo under reasonable assumptions 
for the structure of the intracluster magnetic field. The mild 
decrease of the magnetic field towards the Coma center, which 
reproduces the radial trend of the observed RM distribution in 
Coma, could be better tested with a larger dataset of Faraday 
rotation measures of background radio sources obtainable with 
the next generation sensitive radio telescopes (LOFAR, SKA), 
and with the help of numerical MHD simulations. Radio data 
are the main constraint, so far, to WIMP models. 

The ICS emission produced by the same secondary elec- 
trons is able, in principle, to reproduce both the spectrum and 
the spatial distribution of the EUV emission observed in Coma, 
provided that a quite small average magnetic field B ~ 0.15 
/iG is assumed. Such low value of the B field is also able to 
make the radio data and the hard X-ray data of Coma consis- 
tent within a Synchrotron/IC model for their origins. However, 
such low magnetic field values in Coma produce an unaccept- 
ably large gamma-ray flux, which exceeds the EGRET upper 
limit. The gamma-ray constraints are thus the most stringent 
ones for the analysis of the astrophysical features of DM anni- 
hilations. 

In conclusion, the viable models of WIMP annihilation 
which are consistent with the available data for Coma yield 
a nice fit to the radio data but produce relatively low inten- 
sity emission at EUV, X-ray and gamma-ray frequencies. The 
hadronic gamma-ray emission could be, nonetheless, detected 
by the GLAST-LAT detector. These models also produce neg- 
ligible heating rates for the kind of non-singular halo profile 
we worked out in this paper. It is interesting that the best-fit 
(66) WIMP model with M x = 40 GeV predicts a detectable 
SZ effect (with a peculiar spectrum very different from that 
of the thermal SZ effect) at the level of ~ 40 to 10 /jK in 
the frequency range ~ 10 — 200 GHz, which could be ob- 
servable with the next generation high-sensitivity bolometric 
arrays, space and balloon-borne microwave experiments, like 
PLANCK, OLIMPO, APEX, ALMA. 

The observational "panorama" offered by the next com- 
ing radio, SZ, and gamma-ray astronomical experiments might 
produce further constraints on the viable SUSY model for 
Coma and for other large-scale cosmic structures. Direct DM 
detection experiments have already explored large regions of 
the most optimistic SUSY models, and the planned increase 
in sensitivity of the next-generation experiments will probably 
be able to explore even the core of the SUSY models. In this 
context, we have shown that indirect DM detection proves to 
be not only complementary, but also hardly competitive, es- 
pecially when a full mult i- frequency approach is chosen. When 



combined with future accelerator results, such multi-frequency 
astrophysical searches might greatly help us to unveil the elu- 
sive nature of dark matter. 
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Appendix A: A solution to the diffusion equation 

To understand quantitatively the role of the various popula- 
tions of secondary particles emitting in the Coma cluster, we 
have to describe in details their transport, diffusion and energy 
loss. We consider the following diffusion equation (i.e. neglect- 
ing convection and re-acceleration effects): 



d dn e 
dt~d~E 



D(E,x)V 



+Q e (E,x) 



dn e 
IE 



+ 



d_ 

dE 



b(E,x 



dn e 

He 



(A.l) 



We search for an analytic solution of the diffusion equation in 
the case of diffusion coefficient and energy loss term that do 
not depend on the spatial coordinates, i.e. we take: 



D = D{E) 
b = b(E) 



(A.2) 
(A.3) 



and we implement a slight variant of the method introduced 
in |Baltz fc Edsjo 1998| and IBaltz fc Wai 20041 
Let us define the variable u as: 



y ' dE du 
which yields 



dE' 



(A.4) 



(A.5) 



Then, it follows that b(E) — E/ti oss in terms of the time scale 
ti oss for the energy loss of the relativistic particles, which, for 
E mBX = oo, gives u = t. 
The diffusion equation can be rewritten as 



^=b(E)Q e (E,x) 
du 



(A.6) 



We search for the Green function G of the operator on the left- 
hand-side. Consider the equation for its 4-dimensional Fourier 
transform (t — » u, x — > k): 



-iuj + D(E)k 2 - -jL 
ou. 

which has the solution 
1 
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Transforming back from the Fourier space we find: 

\x — x 

exp 



1 
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where we denned dv = D(u)du, i.e. v — duD(u). The 

suffix 'free' refers to the fact that there are no boundary con- 
ditions yet. These are implemented with the image charges 
method. To apply this technique to galaxy clusters, we can 
consider the approximation of spherical symmetry with Green 
function vanishing at the radius rn- Introducing the set of im- 
age charges (r n ,9 n ,(p n ) = ((— l) n r + 2nrh,d, (j)), one can verify 
that 



the discussion. To do that, we need to specify D(E) and b(E). 
For the diffusion coefficient we assume the form: 



D(E) = Do 



d 2/3 

B 1 ' 3 V 1 GeV 



(t 



1/3 
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(A.15) 



G(r,Y) = 



E 



(-l)"G frcc (r„,Y) 



(A.10) 



fulfills such boundary condition (here Y labels the other vari- 
ables in the Green function) . Moreover, we choose the reference 
frame in such way that we look at the signal along the z polar 
axis (cos 9 = 1) so that \x' — x n \' 2 — (r 1 ) 2 + r 2 — 2 cos 0'r'r n . If 
the source function does not depend on 9' and <j>' , the integral 
on these two variables can be performed explicitly and we find 



(IColafra nccsco & Bla si 19981 IBlasi fc Cola francesc o 1999t 
where ds is the minimum scale of uniformity of the magnetic 
field in kpc (throughout the paper we assume ds — 20 for 
Coma), _B M is the average magnetic field in fiG units, and Do 
some constant that we estimate as Do = 3.1 x 10 28 cm 2 s _1 . 

The energy loss term is the sum of effects due to 
Inverse Compton, synchrotron radiation, Coulomb losses and 
Bremsstrahlung: 

b(E) = b IC (E) + b syn (E) + b C oui(E) + b brern (E) 
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(A.16) 



with t' = t — [u — u') (or no time dependence for stationary 
source). Note that E' > E (energy is lost) and hence u < u, 
v' < v and t' < t. 

A.l. Stationary limit and role of spatial diffusion in 
Coma 

In the limit of time-independence of the source and elec- 
tron number density that has already reached equilibrium, 
Eq. JA.llll takes the form: 



Here n is the mean number density of thermal electrons in 
cm -3 (see Eq. HSU , the average over space gives about n ~ 
3 10~ 3 ), 7 = E/m e and we find b° IC ~ 0.25, b° syn ~ 0.0254, 
uq ou1 ~ 6.13 and b® rt , m ~ 1.51, all in units of 10 -16 GeVs -1 . 
For GeV electrons and positrons the Inverse Compton and syn- 
chrotron terms dominate (see Fig lA.3t . 

To get a feeling about what is the electron/positron en- 
ergy range which will be of interest when considering the radio 
emissivity, we can resort to the "monochromatic" approxima- 
tion, with relativistic particles of a given energy E radiating 
at a single frequency, namely the peak frequency: 
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In the limit in which electrons and positrons lose energy on a 
timescale much shorter than the timescale for spatial diffusion, 
i. e. if the first term on the r.h.s. of Eq. I|A,10 can be neglected, 
the expression for equilibrium number density becomes: 



(dE) ns / T,E) b(E) 



dE' Qe(r,E') 



(A.14) 



This is analogous to the form in Eq. 1A.12L except for the 
factor G (r,v — v') in the integrand: it follows that the latter is 
the Green function term which we need to study to understand 
whether spatial diffusion is important or not. 

Since we have encoded the dependence on the energy loss 
term and the diffusion coefficient in the definition of the vari- 
able v, we preliminarily study what range of Av is relevant in 



Since radio data on Coma extend down to about 30MHz, for 
magnetic fields not much larger than 10 /iG, this translates 
into radiating particles with energies larger than about 1 GeV. 

As a sample case, in Fig. EH we consider a WIMP mass 
M x = 100 GeV and sketch the mapping between the energy 
E' € (E,M X ), with E some reference energy after diffusion, 
and the square root of Av = v — v' = v(E) — v(E'), for a few 
values of E and of the mean magnetic field S M . We find as 
v^" largest value (Av) 1 '" 2 ~ 35 kpc, corresponding to E = 1 GeV 
and _B M = 1 [iG; the maximum value of Av diminishes rapidly 
2 (r') when increasing E or _B M . 

On the right-hand side of Fig. IA.1I we plot G as a function 
of (Av) 1 / 2 for a few values of the radial coordinate r and in 
case the DM halo of Coma is described by a N04 profile. In the 
very central part of the halo, there are significant departures 
of the value of G from unity, on scales (Au) 1//2 at which, for 
the given radius r, the mean squared value of the DM profile 
is significantly different from the square of the value of the 
profile at r. Note, however, that this effect is confined in the 
innermost region of the cluster, corresponding to an angular 
size of ~ 1 — 2 arcmin. We then expect that taking into account 
spatial diffusion will modify only slightly the predictions for the 
radio surface brightness distribution from moderate to large 
radial distances in Coma. 

In Fig. rO we plot G as a function of (Av) 1 / 2 for the same 
values of the radial coordinate r as in Fig. rO but now m case 
the DM halo of Coma is described by a Burkert profile. It is 
clear that departures from unity are essentially negligible even 
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Fig. A.l. Left: The figure shows the distance (Aw) 1 / 2 which, on average, a electron covers while losing energy from 
its energy at emission E' and the energy when it interacts E, for a few values of E: 30 GeV, 10 GeV, 5 GeV and 
1 GeV, and for a few values of the magnetic field (in /iG); we are focusing on a WIMP of mass 100 GeV, hence cutting 
E' < 100 GeV. Right: Green function G as a function of (Au) 1 / 2 , for a few values of the radial coordinate r (in kpc) 
and in case the DM halo of Coma is described by a N04 profile with M v - lr = 0.9 10 15 M Q h~ 1 and c V i r — 10. 



^ 10 

< 




10 



10 



Burkert profile 



200 / 500 



r=l/5/10 



10 15 20 25 30 35 40 

(Av^fkpc] 



10 

< 



10 



10 



: M05 profile 


1 1 1 1 1 1 1 1 I 1 1 l ; 

100 


irv~\ 50^-— 


— ^^200 




500 : 


-\ \ 




1 \ X 




■ \r=l \ 

i , A, , , , 1 , , , , 1 , i , , 1 , , 





10 15 20 25 30 35 



40 



(Av) 1/2 [kpc] 



Fig. A. 2. Green function G as a function of (Au) 1 / 2 , for a few values of the radial coordinate r (in kpc) and in case 
the DM halo of Coma is described by a cored Burkert profile with M V i r — 0.9 10 15 M Q and c„j r = 10 and for the 
Dicmand et al. profile with the same parameters. 



in the inner portion of the halo and, hence, spatial diffusion 
can be safely neglected for all practical purposes in this case. 

To get a more physical insight on the reason why spatial 
diffusion can be neglected, it is useful to consider the follow- 
ing qualitative solution (see, e.g.. ICol afranccsco 2005J for the 
average electron density 



dn e (E,r) 
dE 



[Q e {E,r)Ti oas ] X 



Vs 



TD 



V s +Vo TD + Tl oss 



(A.18) 



which resumes the relevant aspects of the transport equation 
fEa.rO). Here, V s oc R 3 h and V„ oc A 3 (£) are the volumes oc- 
cupied by the DM source and the one occupied by the diffusing 
electrons which travel a distance A(_E) w [D(E) ■ r; oss (ii/)] 1 '' 2 
before loosing much of their initial energy. The relevant time 
scales in Eq. (|A. ip are the diffusion time-scale, td ~ Rf i /D(E), 
and the energy loss time-scale t; oss = E/b e (E), where D is 
again the diffusion coefficient for which we can assume the 
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Fig. A. 3. A comparison among the time scales for the 
energy losses due to various mechanisms (as labeled in the 
figure) and the time scale for diffusion (black solid curve) 
in a cluster of size = 1 Mpc. A uniform magnetic field 
of value B = f /iG and a thermal gas density n — f .3 • f -3 
cm -3 bave been assumed in the computations. 

generic scaling D(E) = D (E / E )~< B" 7 , and b(E) the elec- 
tron energy loss per unit time at energy E. 

For E > E* = {b Eo/R 2 h b B^) ini -^ (for simplic- 
ity we have kept leading terms only, implementing b(E) ~ 
b a (B,j,)(E/GeV) 2 + b Co ui), the condition td > n osa (and con- 
sistently \(E) < Rh) holds, the diffusion is not relevant and 
the solution of Eq. I|A.10 is dn e /dE ~ Q e (E,r)Tioss and shows 
an energy spectrum ~ Q(E) ■ E~ x . This situation (\{E) < R h , 
td > rioss) applies to the regime of galaxy clusters which we 
discuss here for the specific case of Coma, as one can see from 
Fig. ITOl 

For E < E t , the condition td < ti 033 (and consistently 
X(E) > Rh) holds, the diffusion is relevant and the solution 
of Eq. 1A.11 is dn e /dE ~ [Q e {E, t)td] X (V s /V ) and shows an 
energy spectrum ~ Q(E) ■ i5( 2_5 T)/ 2 which is flatter or equal 
to the previous case for reasonable values 7 = 1/3 — 1. This 
last situation {X(E) > Rh, td < n oss ) applies to the regime of 
dwarf galaxies and we will discuss this case more specifically 
elsewhere ( Colafranc esco, Profumo fc Ullio 2 006 ) . 

Fig lA.4l shows the energy shape of the electron equilib- 
rium spectra derived in our approach for a (bb) model with 
M x = 40 GeV and for a W + W~ model with M x = 81 GeV. 
The astrophysical predictions of these two models will be ex- 
tensively discussed in the following. We notice that the energy 
losses in the diffusion equation erase almost completely the 
details of the electron source spectra (see Fig. 0. The equi- 
librium spectra are generally characterized by three different 
regions: i) a low-energy plateau at E i$ 0.1 GeV with a con- 
stant value of dn e /dE which remains almost constant down to 
the electron rest-mass energy; ii) an almost power-law branch 
at 0.1M X ;$ E 0.5M X which is steeper in the softer bb an- 
nihilation final state with respect to the hard spectrum due 
to a W + W~ channel; and iii) a sharp cut-off at the energy 
corresponding to the neutralino mass which marks the natural 
maximum energy of the secondary electron spectra. We will 
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Fig. A. 4. The electron equilibrium spectra calculated at 
the center of Coma as obtained for a soft spectrum due to 
a bb annihilation final state (solid line, model with M x = 
40 GeV) and of a hard spectrum due to a W + W~ channel 
(dashed line, model with M x — 81 GeV). 

show in the next Sect^Jhow these three branches of the elec- 
tron equilibrium spectra will provide observable features in the 
multi-frequency spectrum of Coma and can, consequently, be 
used to constrain the neutralino model. 
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